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]);}]
DSolve
instead ofNDSolve
works for me after cleaning up your ode. But I don't see that you are saving any results in yourFor
loop. $\endgroup$