Numerical continuation methods for bypassing a singularity when integrating an ODE

Question

Here is the ODE I want to numerically integrate,

odey=-l (1 + l) R[y] + (k - y) (-2 Derivative[1][R][y] + (k - y) (R^\[Prime]\[Prime])[y])

If we rearrange it in the standard form,

$$R''(y)-\frac{2}{k-y}R'(y)-\frac{l(l+1)}{(k-y)^{2}}R(y)=0$$

we see that it has a regular singular point at $y=k$ where $k<0$. I also attempted to solve this using DSolve, and it gave me a solution that conforms with the form I expect, that is $R(y)=c_{1}(y-k)^{a}+c_{2}(y-k)^{b}$

DSolve[odey == 0, R[y], y]

Now, I want $R(y)$ to vanish for very large $y$ (boundary condition). So the relevant solution (and its derivative) is

R[y_, l_] := 1/(y - k)^((1 + Sqrt[1 + 4 l (l + 1)])/2)
dR[y_, l_] := -(1/2) (1 + Sqrt[(1 + 2 l)^2]) (-k + y)^(-(3/2) - 1/2 Sqrt[(1 + 2 l)^2])

K[q_] := (Sqrt[\[Pi]] Gamma[1/(q - 1)])/((1 - q) Gamma[1/2 ((q + 1)/(q - 1))])
k = K[-2];

where I just set the constant to 1. Now, I get my initial condition from R itself and then numerically solve the ODE across the region $yM$ to $y0$,

el=0;
rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 200, MaxSteps -> 10000};
yM = - 10^3;
y0 = 1; 
R0m = Rationalize[R[yM, el], rat];
dR0m = Rationalize[dR[yM, el], rat];
BCm = {R[yM] == R0m, R'[yM] == dR0m};
EQm = {(odey /. l -> el) == 0};
Rsolm = R /. First@NDSolve[Union[EQm, BCm], {R}, {y, yM, y0}, rules, Method -> "StiffnessSwitching"];

With the range of integration of the ODE, the NDSolve surely fails since it passes through the singular point at $y=k$. I want to solve the ODE for all $l=(0,30)$.

My question is that, is there a convenient numerical continuation method to bypass the singularity so that NDSolve does not fail?

Answer 1

We can take the approach in Extending NDSolve beyond a singularity and bump it up to order 2. The approach is to solve the ODE in projective space, with $R(z) = [p(x):q(z)]$ in projective coordinates. The real affine solution is then given by $R(z) = p(z)/q(z)$. Since we're adding another variable, we have to add an equation here and there. They are all based on some normalization of the projective coordinates. In this case, we chose $$p(z)^2 + q(z)^2 = \text{constant}\,.$$ By differentiating this relation, we can obtain the extra equations that we need.

projODE = {
    EQm[[1, 1]] /. R -> (p[#]/q[#] &) // Together // Numerator,
    D[p[y]^2 + q[y]^2, {y, 2}]} == 0 // Thread;
projICS = BCm /. {
    R[y0_] == r0_ :> {p[y0], q[y0]} == {Numerator[r0], Denominator[r0]},
    R'[y0_] == rp0_ :> 
     ({D[p[y]/q[y], y] == rp0, D[p[y]^2 + q[y]^2, y] == 0} /. y -> y0)};

projSOL = NDSolve[{projODE, projICS}, {p, q}, {y, yM, y0},
   rules, Method -> "StiffnessSwitching"];

The solution p[y]/q[y] tracks Rsolm[y] up to the singularity, and then continues on past it.

Plot[{Rsolm[y], p[y]/q[y] /. First@projSOL}, {y, yM, y0}, 
 PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[2]},
 PlotLegends -> {R, p/q}]

Plot[{Rsolm[y], p[y]/q[y] /. First@projSOL}, {y, -1, y0}, 
 PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[2]}, 
 Exclusions -> (q[y] == 0 /. First@projSOL), PlotLegends -> {R, p/q}]