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

Question

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:

NDSolve`ProcessSolutions::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.

Answer 1

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.