# 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.