@Preben Alsholm

Perhaps an even more important lesson is first check the dimension of an ODE.

@Preben Alsholm 

Dear Preben, I think I made a major mistake. In my DE is missing a square.

ode  := diff(1/r(phi),phi$2) + 1/r(phi) = G*M/h^2 + 3*G*M/(c^2*r(phi)^2);

What a silly mistake. When I run the program now, the constant is indeed appr. 0.33625

Sorry for the confusion

Last part is becoming a mess now. Filling in the data, the output is - at least for me - unreadable. 

@John Fredsted 

The textbook is using M = 1,99*10^30. All values for G, c and are correct (according to the textbook.)

But I think you are right, the Mass of Mercury is indeed 3.3x10^23 Kg

@John Fredsted 

Thanks, I immediatly downloaded the manual (quite a large manual!!!). In earlier days I've worked with mathcad. But this is slightly different!

Last question of this excercise. It is not said use Maple or Matlab, but can it be done with Maple as well?

The orbital velocity of the planets in our solar system is such that relativistic effects are largely negligible, with the exception of the planet Mercury. In SI units, the data for the orbit of Mercury are M = 1.99 × 10^30, G = 6.67 × 10^−11, h = 2.7 × 10^15, c = 3.0 × 10^8, and the initial conditions can be taken to be r(0) = 4.6 × 10^10, r′(0) = 0.

Show that the precession of the perihelion of Mercury, in this
model, is approximately 0.506 × 10^−6 radians per revolution, and
that since Mercury orbits the Sun every 88 days, this equates to
43 seconds of arc per terrestrial century. 

@John Fredsted 

@Preben Alsholm

Both thanks for the lecture in Maple. I am beginning to like it. I just did some copy paste of your solutions. Are there any good courses in Maple? (Online)?

@Preben Alsholm 

I tried a little myself. Can this be correct?

restart:
params:={M = 1,G = 1,h = 1, c=8};
            
ode  := diff(1/r(phi),phi$2) + 1/r(phi) - G*M/h^2 -3*G*M/((c^2)*r(phi));
                  
             
       
res := dsolve({eval(ode,params),r(0)= 2/3,D(r)(0) = 0},numeric);
               
plots:-odeplot(res,[r(phi)*cos(phi),r(phi)*sin(phi)],0..10*Pi,scaling=constrained,labels=[x,y]);

@Preben Alsholm 

To both of you: thanks. Physically I understand what is going on, but I wouldn't have done it my self by maple. Perhaps you can help me with the second part of Einstein's theory, using Maple17.

d²/d(φ)²(1/r(φ)) + 1/r(φ) = GM/h² + 3GM/(c²*r(φ))

Where c is the speed of light in vacuum. 

Try to use Maple to solve this differential equation numerically, taking M = 1, G = 1, c = 8, h = 1 with initial conditions r(0) = 2/3,

r′(0) = 0. Using polar coordinates, create a plot of the orbit
(r(φ),φ) for 0 ≤ φ < 10π.
The path should resemble an ellipse that rotates about one focus,
in this case the origin. The perihelion of the orbit is defined as the
distance from the point of rotation (in this case the origin) to the
closest point on the path. Here the perihelion rotates about the
origin, so an orbital motion such as this is said to have a precessing perihelion. 

By finding the angular positions φ1, φ2, φ3, φ4 of the first four perihelia, show that the perihelion precesses by approximately a constant amount per revolution δφ = φi+1 − φi − 2π ≃ 0.336 25 radians. 

Thx in advance for the help!