# NDSolve and interpolating function

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[][\[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 of NDSolve works for me after cleaning up your ode. But I don't see that you are saving any results in your For loop. Feb 1, 2019 at 1:21
• The for loop should serve as interpolating function I think? The ode that enters in the for loop is $\textbf{odey}$, which is the nondimensionalized form of the ode. If I run the for loop, it should not return any output. But in my case, there are many errors. What did you do to fix it up? I should use NDSolve since I want to get the numerical solutions. Feb 1, 2019 at 1:42
• I assume you really mean R'[p] and R''[p] in your ode. I did not understand Derivative[][p] and MMA didn't like it either. You can plug numbers in results returned by DSolve. Feb 1, 2019 at 1:53
• can i see your code? Feb 1, 2019 at 1:59
• ah yes, it's the first and second derivative of R. Sorry for that Feb 1, 2019 at 2:03

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][] // 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][] // 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.