Question: How do I compute a Numerical Solution to the Orbital Equation?

I am trying to use Maple to compute the numerical solution to the Newton/Einstein orbital equations and then to plot r(phi) against phi in polar coordinates (obviously without the relativistic correction this should be a flawless ellipse).

I have seen that problems trying to solve ODEs/PDEs are fairly common with people who are new to Maple.  The problem I have is that when I define the ODE with the second derivative, it evaluates the derivative so when I try to use dsolve I get the error message that it was 'expecting' ode to be a differential equation.

ode := diff(1/r(phi), [phi$2]) + 1/r(phi)=GM/h^2;     [(where constants are defined such that GM/h^2 = 1)]
dr :=diff(r(phi),phi);

ics := r(0)=2/3, dr(0)=0;


Also, I am not sure how to define one of the conditions when it is a derivative such that it can be included with sys, although I hade made an attempt.  The other examples I have seen just have functions defined at a particular value, rather than a derivative of the function.  When I try to use dsolve (regardless of the other problem), I get the message 'Error, (in dsolve) not a system with respect to the unknowns {r(0), (diff(r(phi), phi))(0)}'.  Once I have the solution I would like to plot r(phi) against phi for 0 < phi < 2pi to observe the full orbit and then do the same for 0 < phi < 10pi once I have repeated the numerical solution for the orbital equation with relativistic correction so that I can observe that the precession of the perihelion between each revolution is approximately constant. 

I am also not sure how to find the angular positions of the four perihelia after the first revolution, is there a command that will enable me to determine the distance from the origin of the orbit to the closest point on the orbital path for all of the revolutions over the range 0 < phi < 10pi, as I will then be able to make a comparison with each one to check that the precession is approximately constant.

Please Wait...