# NDSolve::ndsz: step size is effectively zero; singularity or stiff system suspected + other warnings for system of differential equations

I'm trying to solve a set of differential equations numerically to get a 3D plot, but I am getting multiple different warnings and errors. First of all, here is the code:

ClearAll["Global*"]
tdot[\[Lambda]_, q_] := Simplify[1/(1 - 1/r)] /. r -> r[\[Lambda]]
rdot[\[Lambda]_, q_] := Simplify[(1/r)*Sqrt[r^2 - q^2*(1 - 1/r)]] /. r -> r[\[Lambda]]
\[Phi]dot[\[Lambda]_, q_] := Simplify[-(q/r^2)] /. r -> r[\[Lambda]]
tasym[\[Lambda]_, q_] := \[Lambda] + Log[\[Lambda]] - 1/\[Lambda] + (q^2 - 2)/(4*\[Lambda]^2) + (3*q^2 - 4)/(12*\[Lambda]^3) + (-3*q^4 + 8*q^2 - 8)/(32*\[Lambda]^4) +
(-9*q^4 + 20*q^2 - 16)/(80*\[Lambda]^5)
rasym[\[Lambda]_, q_] := \[Lambda] + q^2/(2*\[Lambda]) - q^2/(4*\[Lambda]^2) - q^4/(8*\[Lambda]^3) + (3*q^4)/(16*\[Lambda]^4) + (q^6/16 - q^4/20)/\[Lambda]^5
\[Phi]asym[\[Lambda]_, q_] := q/\[Lambda] - q^3/(3*\[Lambda]^3) + q^3/(8*\[Lambda]^4) + q^5/(5*\[Lambda]^5)
\[Lambda]min = 0;
\[Lambda]max = 20;
\[Lambda]inf = 999;
qlist = Array[N[#1] & , 100, {-20, 20}];
maxder = 999999;
eq[q_] := {t[\[Lambda]], r[\[Lambda]], \[Phi][\[Lambda]]} /. NDSolve[{Derivative[1][t][\[Lambda]] == tdot[\[Lambda], q], Derivative[1][r][\[Lambda]] == rdot[\[Lambda], q],
Derivative[1][\[Phi]][\[Lambda]] == \[Phi]dot[\[Lambda], q], t[\[Lambda]inf] == tasym[\[Lambda]inf, q], r[\[Lambda]inf] == rasym[\[Lambda]inf, q],
\[Phi][\[Lambda]inf] == \[Phi]asym[\[Lambda]inf, q], WhenEvent[{Abs[Derivative[1][t][\[Lambda]]] > maxder ||
Abs[Derivative[1][r][\[Lambda]]] > maxder || Abs[Derivative[1][\[Phi]][\[Lambda]]] > maxder}, "StopIntegration"]},
{t[\[Lambda]], r[\[Lambda]], \[Phi][\[Lambda]]}, {\[Lambda], \[Lambda]min, \[Lambda]inf},
{"ExtrapolationHandler" -> {Indeterminate & , "WarningMessage" -> False}}][[1]]
eqlist = (eq[#1] & ) /@ qlist;
tlist = eqlist[[All,1]];
rlist = eqlist[[All,2]];
\[Phi]list = eqlist[[All,3]];
surface = MapThread[{#2*Sin[#3], #2*Cos[#3], If[Abs[#3] < Pi, #1, Indeterminate]} & , {tlist, rlist, \[Phi]list}];
Show[ParametricPlot3D[surface, {\[Lambda], \[Lambda]min, \[Lambda]max}, PlotRange -> {{-20, 20}, {-30, 20}, {-30, 20}}], ImageSize -> Large]


Running this code, with the parameters I chose after playing around with them, gives no warning message, and results in what i'm trying to get, but only half of it. For certain values of the parameters $$\lambda_{inf}$$, $$maxder$$, the warning message in the title of the question appears, however after reading other questions regarding this issue I don't think this is too much of a problem.

The main problems start when i set $$\lambda_{min}=-30$$ which would give me the other half of the plot. the first warning message that appears is

NDSolve::mxst: Maximum number of 129336 steps reached at the point [Lambda] == -0.53469.

same thing for other values of $$\lambda$$. At first I tried overcoming this by increasing MaxSteps, however this didn't work for me as Mathematica would just use up all my RAM and my computer would stop working.

To investigate I tried to solve only for the negative values, I set $$\lambda_{min}=-30$$ and $$\lambda_{max}=-5$$, and set NDSolve to solve only between $$\lambda_{min}$$ and $$\lambda_{max}$$. What happens with these settings is the following warning message:

NDSolveProcessSolutions::nodata: No solution data was computed between [Lambda] == -30. and [Lambda] == -5..

Which is weird since I should have a solution, unless I made some dumb mistake. Reading other questions, I saw this could maybe be a case of backslide where since the solution doesn't adhere to the new standards of NDSolve it doesn't give any solution. However this is just speculation.

I also tried adding some methods to NDSolve with no results, since I don't know much about them. What I hope is that by tweaking some NDSolve parameters or using some methods I can manage to get a result, so any suggestion in this direction is very welcome.

Rationalize[qlist] and ParametricNDSolveValue evaluates without error

 tr\[Phi] =ParametricNDSolveValue[{Derivative[1][t][\[Lambda]] ==tdot[\[Lambda], q],Derivative[1][r][\[Lambda]] == rdot[\[Lambda], q],Derivative[1][\[Phi]][\[Lambda]] == \[Phi]dot[\[Lambda], q],t[\[Lambda]inf] == tasym[\[Lambda]inf, q],r[\[Lambda]inf] ==rasym[\[Lambda]inf, q], \[Phi][\[Lambda]inf] == \[Phi]asym[\[Lambda]inf, q] ,WhenEvent[{Abs[Derivative[1][t][\[Lambda]]] > maxder ||Abs[Derivative[1][r][\[Lambda]]] > maxder ||Abs[Derivative[1][\[Phi]][\[Lambda]]] > maxder}, "StopIntegration"] }
, {t ,r , \[Phi] }, {\[Lambda], \[Lambda]min, \[Lambda]inf}  , {q}, \{"ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}}]

Plot[Table[Through[ tr\[Phi][q][\[Lambda]]],{q,Rationalize[qlist]}]  , {\[Lambda], \[Lambda]min, \[Lambda]inf},Evaluated -> True]


WhenEvent seems to be unnecessary.

• This does not seem to work for me for the same reason my code doesn't work, when i put $\lambda_{min}=-30$ the domain of the solution is $[-0.33,\lambda_{inf}]$ but I'm expecting a solution also for lower values of $\lambda$. Jun 17, 2021 at 12:17
• tr\[Phi][-30 ] gives a solution in the range 0<\[Lambda]<999 Jun 17, 2021 at 12:26
• Yes, but i'm also looking for $-30<\lambda<0$, not q=-30. That is my main problem, my original code already gives a good plot for $0<\lambda<999$. But Mathematica says there is no solution for that domain, or reaches MaxSteps and stops. Maybe i could force NDSolve to be less precise somehow? Jun 17, 2021 at 12:34
• Try Method->"StiffnessSwitching" Jun 17, 2021 at 12:51
• Thank you, combining this with your code, and some tweaking on precision and accuracy seems to give a result. I'm going to try and get a plot out of it now Jun 17, 2021 at 13:00