NDSolve Maxsteps Reached Error

Question

I'm trying to solve a system of Equations. The code is commented below, I am mainly conserved with the last few lines, specifically the NDSolve line. No matter what I do to the Maxsteps option, the limit is reached at the same point t=16.486...

I'm not sure why NDSolve stops at this point. Any help would be greatly appreciated!

Edit: I've attempted to use AccuracyGoal and PrecisionGoal, they've stopped the error of maxsteps, but when I plot r(t) it only plots up until the time limit above.

(*all quantities are in their natural units, per unit kg*)
(* angular momentum of star in J (conserved i.e constant) *)
L = 3.086*10^25;
(* V is a constant in km/s *)
V = 200*10^3;
(* initial radial velocity (in polar coordinates) *)
vr = 50*10^3;
(* initial tangential velocity (in polar coordinates) *)
vt = 10^5;
(* initial radius *)
r0 = 10*3.086*10^19;
(* seconds in a mega year (units of time I want to plot in *)
megayear = 3.154*10^7*10^6;
(* energy of the star system calculated from initial conditions \
(conserved i.e constant energy), the V^2 Log(r0) term is the \
spherical potential *)
energy = vr^2/2 + L^2/(2*r0^2) + V^2*Log[r0];
r =.
(* solve for the max and min radii of the orbit *)
Solve[L^2/(2 r^2) + V^2 Log[r] == energy, r]

(* Problem = If I increase the range over t (say, 100 instead of 10) \
I get a maxstep error *)
s = NDSolve[{r'[t] == 
    Sqrt[2*(energy - V^2*Log[r[t]] - L^2/(2*r[t]^2))*
      megayear^2], \[Theta]'[t] == L/r[t]^2*megayear, \[Theta][0] == 
    vt/r0, r[0] == r0}, {r, \[Theta]}, {t, 0, 100}, MaxSteps -> 100]

Answer 1

Sign r'[t] changes when reflected from points rm = r /. Solve[L^2/(2 r^2) + V^2 Log[r] == energy-vr^2/2, r]. Therefore, the model should be subject to the sign r'[t]. But it’s better to use a second order equation

(*all quantities are in their natural units,per unit kg*)(*angular \
momentum of star in J (conserved i.e constant)*)L = 3.086*10^25;
(*V is a constant in km/s*)
V = 200*10^3;
(*initial radial velocity (in polar coordinates)*)
vr = 50*10^3;
(*initial tangential velocity (in polar coordinates)*)
vt = 10^5;
(*initial radius*)
r0 = 10*3.086*10^19;
(*seconds in a mega year (units of time I want to plot in*)
megayear = 3.154*10^7*10^6;
(*energy of the star system calculated from initial conditions \
(conserved i.e constant energy),the V^2 Log(r0) term is the spherical \
potential*)
energy = vr^2/2 + L^2/(2*r0^2) + V^2*Log[r0];
rm = r /. 
  Solve[L^2/(2 r^2) + V^2 Log[r] == energy - vr^2/2, 
   r];(*solve for the max and min radii of the orbit*)


s = NDSolve[{r''[t] == (-V^2/r[t] + L^2/r[t]^3)*megayear^2, \[Theta][
     0] == vt/r0, r[0] == r0, r'[0] == vr}, {r, \[Theta]}, {t, 0, 
   400}]
Plot[{r[t] /. s, rm}, {t, 0, 400}, PlotRange -> All, 
 AxesOrigin -> {0, 0}, AxesLabel -> Automatic]