# Numerical continuation methods for bypassing a singularity when integrating an ODE

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?

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}]