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?