NDSolve and interpolating function

Question

What could possibly went wrong in my code? Basically, I am solving the differential equation $\textbf{ode}$ using $\textbf{NDSolve}$. But mathematica says, NDSolve::mxst: Maximum number of 10000 steps reached at the point y == -0.431184926619225550730163530627, and InterpolatingFunction::dmval: Input value {1/1000} lies outside the range of data in the interpolating function. Extrapolation will be used.

What might cause these errors and how to address these? Also, there is an infinite expression $\frac{1}{0}$ encountered in my code. Is there a better way to execute my code?

 rl[\[Rho]_, b0_, q_] := \[Rho] (1 - (Sqrt[\[Pi]]Gamma[1/(q - 1)])/((1 - q) Gamma[1/2 ((q + 1)/(q - 1))])b0/\[Rho])
 r1 = FullSimplify[rl[\[Rho], a, -2]];
 f0 = (1 - a^3/r1^3) /. \[Rho] -> a y0;

 ode = -l (1 + l) R[\[Rho]] + 2 (\[Rho] + (a Sqrt[\[Pi]] Gamma[2/3])/Gamma[1/6]) Derivative[1][][\[Rho]] + (\[Rho] + (a Sqrt[\[Pi]] Gamma[2/3])/Gamma[1/6])^2 (R^\[Prime]\[Prime])[\[Rho]]
 odey = ode /. R -> (R[#/a] &) /. \[Rho] -> a y // Simplify

 R[y_, l_] := (1/y)^(1 + l) - ((2 + 3 l + l^2) \[Pi] (1/y)^(3 + l) Gamma[2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) + ((24 + 50 l + 35 l^2 + 10 l^3 + l^4) \[Pi]^2 (1/y)^(5 + l) Gamma[2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4); 
 dR[y_, l_] := -(1 + l) (1/y)^(2 + l) + ((6 + 11 l + 6 l^2 + l^3) \[Pi] (1/y)^(4 + l)Gamma[2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) - ((120 + 274 l + 225 l^2 + 85 l^3 + 15 l^4 + l^5) \[Pi]^2 (1/y)^(6 + l) Gamma[2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4);

 rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 30, MaxSteps -> 10000};
 rat = 10^-30;
 yP = 10^3;
 yM = -10^3;

 a = 1;
 Q = 1;
 f[r_] := 1 - a^3/r^3;
 \[Psi][r_] := 0;
 aRP = -1/(r^2 Sqrt[f[r]]);
 bRP = -1/(2 r^2) (1 + r \[Psi]'[r]);
 dRP = -1/(16 r^2) ((1 + r \[Psi]'[r]) - (1 + r \[Psi]'[r] + 3 r^2 \[Psi]'[r]^2 + r^3 \[Psi]'[r]^3 - 6 r^2 \[Psi]''[r] - 2 r^3 D[\[Psi][r], {r, 3}]) f[r] + (1 + 4 r \[Psi]'[r] + 3 r^2 \[Psi]''[r]) r f'[r] + (1 + r \[Psi]'[r]) r^2 f''[r]);

 y0 = 10^-3;

 For[el = 0, el <= 5, 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)])/(r1 /. \[Rho] -> a y0)^2,rat];

c1f = (rm s)/(drp rm - drm rp);(*jump condition*)
baremode[el] = Sqrt[(2 el + 1)/(4 \[Pi])] (c1f Rsolp'[y0]);}]

Answer 1

Clear["Global`*"]

rl[ρ_, b0_, 
  q_] := ρ (1 - (Sqrt[π] Gamma[1/(q - 1)])/((1 - q) Gamma[
         1/2 ((q + 1)/(q - 1))]) b0/ρ)
r1 = FullSimplify[rl[ρ, a, -2]];
f0 = (1 - a^3/r1^3) /. ρ -> a y0;

ode = -l (1 + l) R[ρ] + 
   2 (ρ + (a Sqrt[π] Gamma[2/3])/
      Gamma[1/6]) R'[ρ] + (ρ + (a Sqrt[π] Gamma[2/3])/
      Gamma[1/6])^2 R''[ρ];

odey = ode /. R -> (R[#/a] &) /. ρ -> a y // Simplify;

R[y_, l_] := (1/y)^(1 + 
      l) - ((2 + 3 l + l^2) π (1/y)^(3 + l) Gamma[2/3]^2)/(2 (3 + 
        2 l) Gamma[1/6]^2) + ((24 + 50 l + 35 l^2 + 10 l^3 + 
        l^4) π^2 (1/y)^(5 + l) Gamma[2/3]^4)/(8 (15 + 16 l + 
        4 l^2) Gamma[1/6]^4);
dR[y_, l_] := -(1 + l) (1/y)^(2 + 
       l) + ((6 + 11 l + 6 l^2 + l^3) π (1/y)^(4 + l) Gamma[
        2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) - ((120 + 274 l + 
        225 l^2 + 85 l^3 + 15 l^4 + l^5) π^2 (1/y)^(6 + l) Gamma[
        2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4);

rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, 
   WorkingPrecision -> 30, MaxSteps -> 10000};
rat = 10^-30;
yP = 10^3;
yM = -10^3;

a = 1;
Q = 1;
f[r_] := 1 - a^3/r^3;
ψ[r_] := 0;
aRP = -1/(r^2 Sqrt[f[r]]);
bRP = -1/(2 r^2) (1 + r ψ'[r]);
dRP = -1/(16 r^2) ((1 + 
       r ψ'[r]) - (1 + r ψ'[r] + 3 r^2 ψ'[r]^2 + 
        r^3 ψ'[r]^3 - 6 r^2 ψ''[r] - 
        2 r^3 D[ψ[r], {r, 3}]) f[
       r] + (1 + 4 r ψ'[r] + 3 r^2 ψ''[r]) r f'[
       r] + (1 + r ψ'[r]) r^2 f''[r]);

y0 = 10^-3;

el = 0

dR0p = Rationalize[dR[yP, el], rat];
R0m = Rationalize[R[yM, el], rat];
dR0m = Rationalize[dR[yM, el], rat];

Rsolp[y_] = 
  R[y] /. DSolve[{(odey /. l -> el) == 0, R[yP] == R0p, 
       R'[yP] == dR0p}, R[y], y][[1]] // Simplify;

Rsolp[y] // N // Simplify

(*(1.00086 - 4.31061*10^-7 y)/(1. y + 0.431185)*)

Rsolm[y_] = 
  R[y] /. DSolve[{(odey /. l -> el) == 0, R[yM] == R0m, 
       R'[yM] == dR0m}, R[y], y][[1]] // Simplify;

Rsolm[y] // N // Simplify

(*(0.999137 - 4.31309*10^-7 y)/(1. y + 0.431185)*)

The you can set a new el and do it again.