The way you are doing it looks pretty good to me. Your residuals are down in the region of 10^-8 to10^-9, and given that your starting point is a numerical solution from odesolve anything better would be meaningless. The only thing I would do differently is that I would set up a vector of times, and then have one vectorized equation in the solve block rather than 12 separate equations.
It will be one chapter of the book "Fine Math with Mathcad". Would you be one coauthor
Thanks for the consideration, but that sounds like it might be a lot of work
I can nor create one orbit as parabola
You can't have an orbit as a parabola, because if it's a parabola it's not an orbit
The parabolic trajectory is the limiting case of the hyperbolic trajectory, when the orbiting body just achieves escape velocity ( i.e. it's the boundary case between an elliptical orbit and a hyperbolic trajectory). In such a case there is a very specific relationship between the orbital velocity and the radial distance. See the attached worksheet
