This notebook plots the orbits of some central force problems.  Here is a general ODE integrating function which calls Mathematica's built-in NDSolve.

\(Orbit[ode_, {x0_, y0_, u0_, v0_}, 
        time_]\  := 
	{x[t], y[t], u[t], 
          v[t]} /. \(NDSolve[
            ode~Join~{x[0] == x0, y[0] == y0, u[0] == u0, 
                v[0] == v0}, {x[t], y[t], u[t], v[t]}, {t, 0, time}, 
            MaxStepFraction -> 1/1000]\)[\([1]\)]\)

To get a problem with a given potential energy U(r) the function F below should be set to compute  -U'[r]/r.  For the Kepler problem U= -1/r so we set F = -1/r^3.  The orbits are the familiar ellipses with focus at the origin.

F[x_, y_]  :=  
      With[{r = Sqrt[x^2 + y^2]}, -1/r^3 ]
\(\(ode\  = \ {\(x'\)[t] == u[t], \(y'\)[t] == 
            v[t], \(u'\)[t] == \ F[x[t], y[t]]*x[t], \(v'\)[t] == \ 
            F[x[t], y[t]]*y[t]};\)
\)
X1  =  Take[Orbit[ode, {1, 0, 0, 0.5}, 10], 2];
X2  =  
        Take[Orbit[ode, {0, 1, 0.5, 0}, 10], 2];
X3  =  
        Take[Orbit[
            ode, {-1/Sqrt[2], -1/Sqrt[2], 
              0.5/Sqrt[2], -0.5/Sqrt[2]}, 10], 
          2];
pl1 = 
        ParametricPlot[{X1, X2, X3}, {t, 0, 10}, AspectRatio -> 1, 
          PlotRange -> {{-1.1, 1.1}, {-1.1, 
                1.1}}];


X1  =  Take[Orbit[ode, {1, 0, 0, 0.25}, 10], 2];
X2  =  
        Take[Orbit[ode, {0, 1, 0.25, 0}, 10], 
          2];
X3  =  
        Take[Orbit[
            ode, {-1/Sqrt[2], -1/Sqrt[2], 
              0.25/Sqrt[2], -0.25/Sqrt[2]}, 10], 
          2];
pl2 = 
        ParametricPlot[{X1, X2, X3}, {t, 0, 10}, AspectRatio -> 1, 
          PlotRange -> {{-1.1, 1.1}, {-1.1, 
                1.1}}];
Show[
        GraphicsArray[{pl1, 
            pl2}]];

It is very atypical to have all orbits periodic like this.  To illustrate we add a small perturbing potential so that now U = -1/r + k/r^2

F[x_, y_]  :=  
      With[{r = Sqrt[x^2 + y^2]}, -1/r^3 - 
          0.1/r^4 ]

ode  =  {x'[t] == u[t], y'[t] == 
            v[t], u'[t] ==  F[x[t], y[t]]*x[t], v'[t] ==  
            F[x[t], y[t]]*y[t]};
X1  =  Take[Orbit[ode, {1.2, 0, 0, 1}, 100], 2];
pl1  =  
        ParametricPlot[{X1}, {t, 0, 100}, AspectRatio -> 1, 
          PlotRange -> All];
X1  =  Take[Orbit[ode, {1, 0, 0, 0.9}, 100], 2];
pl2  =  
        ParametricPlot[{X1}, {t, 0, 100}, AspectRatio -> 1, 
          PlotRange -> All];
X1  =  Take[Orbit[ode, {1, 0, 0, 0.8}, 100], 2];
pl3  = 
        ParametricPlot[{X1}, {t, 0, 100}, AspectRatio -> 1, 
          PlotRange -> All];
Show[GraphicsArray[{pl1, pl2, pl3}]];

You could say that the ellipses are now precesssing.  Note that the rate of precession is not the same for all orbits -- different initial conditions lead to different rates in a given mechanical system.  In other words, there are invariant tori with various frequency ratios in the phase space; of course some of these will be rational and others irrational.

One more example, with U = r^4/r

F[x_, y_]  :=  
      With[{r = 
            Sqrt[x^2 + 
                y^2]}, -r^2 ]

ode  =  {x'[t] == u[t], y'[t] == 
            v[t], u'[t] ==  F[x[t], y[t]]*x[t], v'[t] ==  
            F[x[t], y[t]]*y[t]};
X1  =  Take[Orbit[ode, {1.2, 0, 0, 1}, 100], 2];
pl1  =  
        ParametricPlot[{X1}, {t, 0, 100}, AspectRatio -> 1, 
          PlotRange -> All];
X1  =  Take[Orbit[ode, {1, 0, 0, 0.5}, 100], 2];
pl2  =  
        ParametricPlot[{X1}, {t, 0, 100}, AspectRatio -> 1, 
          PlotRange -> All];
X1  =  Take[Orbit[ode, {1, 0, 0, 0.1}, 100], 2];
pl3  = 
        ParametricPlot[{X1}, {t, 0, 100}, AspectRatio -> 1, 
          PlotRange -> All];
Show[GraphicsArray[{pl1, pl2, pl3}]];