Here are some basic equations
Vn = -Mn/Sqrt[R^2 + z^2 + cn^2]; Vd = -Md/Sqrt[b^2 + R^2 + (a + Sqrt[h^2 + z^2])^2]; Vh = υ0^2/2 Log[R^2 + β*z^2 + ch^2]; Vt = Vn + Vd + Vh; V = Vt + Lz^2/(2*R^2); Vt = V /. {R -> R[t], z -> z[t]}; the values of the parameters
Mn = 250; cn = 0.25; Md = 7000; b = 6; a = 3; h = 0.2; υ0 = 20; β = 1.5; ch = 8.5; Lz = 10; En = 600; the forces along R and z coordinates
FR = -D[V, R]; Fz = -D[V, z]; FRt = FR /. {R -> R[t], z -> z[t]}; Fzt = Fz /. {R -> R[t], z -> z[t]}; the system of the differential equations
DifferentialEquations[R0_, z0_, pR0_, pz0_] := Module[{Deq1, Deq2, Deq3, Deq4}, Deq1 = R'[t] == pR[t]; Deq2 = z'[t] == pz[t]; Deq3 = pR'[t] == FRt; Deq4 = pz'[t] == Fzt; {Deq1, Deq2, Deq3, Deq4, R[0] == R0, z[0] == z0, pR[0] == pR0, pz[0] == pR0} ] the initial conditions
R0 = 5; z0 = 0; pR0 = 0; pz0 = Sqrt[2*(En - V) - pR0^2] /. {R -> R0, z -> z0} tmax = 500; and the integration loop
sdeq = DifferentialEquations[R0, z0, pR0, pz0]; data0 = {}; sol = NDSolve[sdeq, {R[t], z[t], pR[t], pz[t]}, {t, 0, tmax}, MaxSteps -> Infinity, WorkingPrecision -> MachinePrecision, PrecisionGoal -> 12, AccuracyGoal -> 12, Method -> {"EventLocator", "Event" -> z[t], "Direction" -> -1, "EventAction" :> AppendTo[data0, {{R[t], pR[t]}, {R[t], -pR[t]}}], Method -> "Adams"}]; Mathematica does not report any mistakes the integration loop. However no output is produced (empty data0 list) the program consumes all the available memory until it crashes.
Any ideas about what cause this? Also any suggestions for code improvements are more than welcome!
![z[t]](https://i.sstatic.net/ofYbj.png)
R[t]andpR[t]are rapidly oscillating functions, whilez[t]andpz[t]appear to be zero. Note that Fzt is proportional to z, so if z starts at zero it would appear to stay at zero. $\endgroup$data0is not empty? If not, can you please make a quickListPlotof the elements of the list? $\endgroup$data0is empty, as it would be forz[t]identically zero. $\endgroup$