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.