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 I make the transformation $t=1/y$ (where $y=\infty$ corresponds to $t=0$),
odet = odey /. R -> (R[1/#] &) /. y -> 1/t // Simplify
What I do next is to perform Frobenius method (series solution of $R$)
odetSeries = Series[odet /. R -> Function[{t}, t^p Sum[a[i] t^i, {i, 0, 7}]] // Simplify, {t, 0, 4}] // Simplify
where the indicial equation to solve is
Solve[-(l + l^2 + p - p^2) a[0] == 0, p]
I chose $p=1+l$ to satisfy the boundary condition. After which, I obtained the first few coefficients of the expansion $a_{0},a_{1},a_{2}$
odetSeries /. p -> 1 + l /. a[1] -> 0 // Simplify
Solve[(k^2 (2 + 3 l + l^2) a[0] + 2 (3 + 2 l) a[2]) == 0, a[2]] //Flatten[#] & // Simplify
odetSeries /. p -> 1 + l /. a[1] -> 0 /. % /. a[3] -> 0 // Simplify
Solve[(-k^4 (24 + 50 l + 35 l^2 + 10 l^3 + l^4) a[0] + 8 (15 + 16 l + 4 l^2) a[4]) == 0, a[4]] // Flatten[#] & // Simplify
a4[l_] := (k^4 (24 + 50 l + 35 l^2 + 10 l^3 + l^4) a0)/(8 (15 + 16 l + 4 l^2));
a2[l_] := -((k^2 (2 + 3 l + l^2) a0)/(6 + 4 l));
solInfPlus[l_] := t^(1 + l) (a0 + a2[l] t^2 + a4[l] t^4);
testsol = solInfPlus[l] /. t -> 1/y
Collect[testsol, y, Simplify] /. a0 -> 1
D[testsol, y] /. a0 -> 1 // Collect[#, y, Simplify] &
R[y_, l_] := (1/y)^(1 + l) - (k^2 (2 + 3 l + l^2) (1/y)^(3 + l))/(6 + 4 l) + (k^4 (24 + 50 l + 35 l^2 + 10 l^3 + l^4) (1/y)^(5 + l))/(8 (15 + 16 l + 4 l^2))
dR[y_, l_] := -(1 + l) (1/y)^(2 + l) + (k^2 (6 + 11 l + 6 l^2 + l^3) (1/y)^(4 + l))/(6 + 4 l) - (k^4 (120 + 274 l + 225 l^2 + 85 l^3 + 15 l^4 + l^5) (1/y)^(6 + l))/(8 (15 + 16 l + 4 l^2))
Now, I execute numerical calculations.
k=(Sqrt[\[Pi]] Gamma[-(1/3)])/(3 Gamma[1/6]);
rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 30, MaxSteps -> 10000};
rat = 10^-30;
yP = 10^3;
yM = -10^3;
y0 = 10^-2;
For[el = 0, el <= 8, el++, {
R0p = Rationalize[R[yP, el], rat];
dR0p = Rationalize[dR[yP, el], rat];
R0m = Rationalize[R[yM, el], rat];
dR0m = Rationalize[dR[yM, el], rat];
Rsolp = R /. First@NDSolve[{(odey /. l -> el) == 0, R[yP] == R0p, R'[yP] == dR0p}, {R}, {y, yP, y0}, rules];
Rsolm = R /. First@NDSolve[{(odey /. l -> el) == 0, R[yM] == R0m, R'[yM] == dR0m}, {R}, {y, yM, y0}, rules];
rp = Rationalize[ Rsolp[y0], rat];
rm = Rationalize[ Rsolm[y0], rat];
drp = Rationalize[ Rsolp'[y0], rat];
drm = Rationalize[ Rsolm'[y0], rat];
s = Rationalize[(-Q a Sqrt[4 Pi (2 el + 1)])/(a(y0-k))^2,rat];
c1f = (rm s)/(drp rm - drm rp);(*jump condition at y=y0 between Rsolp and Rsolm*)
baremode[el] = Sqrt[(2 el + 1)/(4 \[Pi])] (c1f Rsolp'[y0]);}]
But the solver stops evaluating exactly at $y=k=-0.431188993916544683366628968310$, this is the regular singular point of the ODE above. Specifically, the NDSolve fails in Rsolm since Rsolm covers the region from $yM=-10^{3}$ to $y0=10^{-2}$ in which the singularity $y=k=-0.431188993916544683366628968310$ belongs.
Is there a way to solve this problem? I am thinking of integrating near the singular point, bypass it and then continue integrating. But I do not have an idea on how to execute it properly.
odeyor ` in the standard form`? At first glance, these are different equations. Or did you make a mistake? $\endgroup$