Quasinormal modes with Mathematica

I want to compute the quasi-normal modes for the Kerr-AdS black hole, so basically after decoupling the angular and the radial equation, I have troubles to compute the separation constant and also the quasi-normal modes.

I found on Computation of Kerr QNMs with Leaver’s method a code in mathematica applied to Kerr BH. They follow the Leaver's method, which is based on three-recurrence formulae for the coefficients of the Taylor series solutions. I should get from the angular three-recurrence the separation constant and from the radial part, the frequencies. Unfortunately I can't see it and if someone could help me to understand the code, I will thank you forever.

(We first set the parameters of the Kerr spacetime, the field and the mode we're interested in:)

ar = 0.200;
br = Sqrt[1 - 4*ar^2];
rminus = (1 - br)/2;
rplus = (1 + br)/2; s = -2; m = 2; l = 2;
k1 = 1/2*Abs[m - s];
k2 = 1/2*Abs[m + s];
Alm = l*(l + 1) - s*(s + 1)
(*This is an initialization on Alm*);
winit = 2*(0.37 - 0.08 I);
(*This is an initialization for the frequency*)


(In the following we compute the continued fraction with 130 terms \ (you can do it with more of course) and recursively re-computing the \ angular separation constant. Seven recursions are enough for this \ value of rotation parameter, as seen in the iterated values; See eqs. \ (18)-(27) in [1])

ie = 0;
While[ie < 7, NITMAX = 130; c0 = 1 - s - I*w - 2*I/br*(w/2 - ar*m);
c1 = -4 + 2*Iw(2 + br) + 4*I/br*(w/2 - ar*m);
c2 = s + 3 - 3*I*w - 2*I/br*(w/2 - arm);
c3 = w^2(4 + 2*br - ar^2) - 2*armw - s - 1 + (2 + br)Iw - Alm + (4*w + 2*I)/br*(w/2 - ar*m);
c4 = s + 1 - 2*w^2 - (2*s + 3)Iw - (4*w + 2*I)/br*(w/2 - ar*m);
[Gamma] = Function[n, n^2 + (c2 - 3)*n + c4 - c2 + 2];
[Beta] = Function[n, -2*n^2 + (c1 + 2)*n + c3];
[Alpha] = Function[n, n^2 + (c0 + 1)*n + c0];
Leaver31[w_] := Module[{Rn}, For[{n = NITMAX; Rn = -1.0;}, n > 0, { Rn = [Gamma][n]/([Beta][n] - [Alpha][n]*Rn); n--;} ]; Rn]; Leaver33[w_] := [Beta][0]/[Alpha][0] - Leaver31[w];
wang = FindRoot[Leaver33[w] == 0, {w, winit}][[1]][[2]];
witer[ie] = wang; [Gamma]ang = Function[n, 2*arwang(n + k1 + k2 + s)]; [Beta]ang = Function[n, n*(n - 1) + 2*n*(k1 + k2 + 1 - 2*ar*wang) - (2*ar* wang*(2*k1 + s + 1) - (k1 + k2)*(k1 + k2 + 1)) - (ar^2*wang^2 + s*(s + 1) + Sep)];
[Alpha]ang = Function[n, -2*(n + 1)*(n + 2*k1 + 1)];
Leaver31Ang[Sep_] := Module[{Rn}, For[{n = NITMAX; Rn = -1.0;}, n > 0, { Rn = [Gamma]ang[n]/([Beta]ang[n] - [Alpha]ang[n]*Rn); n--;} ]; Rn]; Leaver33Ang[Sep_] := [Beta]ang[0]/[Alpha]ang[0] - Leaver31Ang[Sep];
Alm = FindRoot[Leaver33Ang[Sep] == 0, {Sep, Alm}][[1]][[2]]; ie = ie + 1]; Table[witer[n], {n, 0, 6}]


First, they define Alm, so I thought they are going to replace this value on the radial three-recurrence to get the frequencies, but I have seen that there is a correction to the quasi normal frequency coming from the angular and also to the separation itself. Does it sound right?

How do I get more frequencies assuming m=0 and l=2 and s=-2? I just get the fundamental one, the other frequencies on the list are the convergence of the ω0 How could I get an analytical expansion for the separation constant using this code?

Thanks!