Classical Kepler Problem

Clear["Global`*"]

<< Graphics`Graphics`

Potential is V(r) = - G m1 m2 / r  [ 1 - G (m1+m2)/(2c^2 r)  + O(1/c^4) ]

v[u_] := -k/r + b/r^2 /. r -> 1/u

D[v[u],u]

u''[phi] + u[phi] == - (m/l^2) D[v[u[phi]],u[phi]] ;

Apart[#] & /@ %

eq = % ;

?/@

k=1; l=1; m=1; b=0;

u0 = 1.6; up0 = 0; n = 10;

eq

sol = NDSolve[ {eq, u[0]==u0, u'[0]==up0},
                u[phi], {phi,0,n 2Pi},
                     MaxSteps -> Infinity] ;

??NDSolve

u1[phi_] = u[phi] /. sol ;

pl1 = PolarPlot[ Evaluate[1/u1[phi]], {phi,0,n 2Pi},
                 PlotStyle->RGBColor[1.0,0.0,0.1] ] ;

?PolarPlot

b = 0.02

eq

sol = NDSolve[ {eq, u[0]==u0, u'[0]==up0},
                u[phi], {phi,0,n 2Pi},
                     MaxSteps -> Infinity] ;

u1[phi_] = u[phi] /. sol ;

pl2 = PolarPlot[ Evaluate[1/u1[phi]], {phi,0,n 2Pi},
                 PlotStyle->RGBColor[0.1,0.1,1.0] ] ;

Show[pl1,pl2] ;