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:

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.

 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]