>E=-0.2;
>l=1;
>
>function phi(r)
$return -1/r;
$endfunction
>
>function dphidr(r)
$return 1/r**2;
$endfunction
>
>function Ueff(r)
$global E,l;
$return phi(r)+l**2/(2*r**2);
$endfunction
>
>r=linspace(0.4,10,96);
>xplot(r,Ueff(r));
>
>function rdot(r)
$global E,l;
$return sqrt(2*E - 2*phi(r) - l**2/r**2);
$endfunction
>
>function rdoubledot(r)
$global E,l;
$return -dphidr(r) + l**2/r**3;
$endfunction
>
>function phidot(r)
$global E,l;
$return l/r**2;
$endfunction
>
>function f(t,x)
$return [x[2] , rdoubledot(x[1]) , phidot(x[1])];
$endfunction
>
>x=[3,rdot(3),0]
            3      0.394405             0 
>t=linspace(0,40,10000);
>sol=runge("f",t,x);
>xplot (t,sol[1]);
>xplot (sol[1]*cos(sol[3]),sol[1]*sin(sol[3]));
>
>
>function phi(r)
$return -1/r**0.95;
$endfunction
>
>function dphidr(r)
$return 0.95/r**1.95;
$endfunction
>
>r=linspace(0.4,10,96);
>xplot(r,Ueff(r));
>
>x=[3,rdot(3),0]
            3      0.439546             0 
>t=linspace(0,80,10000);
>sol=runge("f",t,x);
>xplot (t,sol[1]);
>xplot (sol[1]*cos(sol[3]),sol[1]*sin(sol[3]));
>
>xplot (t,sol[3]);
>2*Pi
      6.28319 
>4*Pi
      12.5664 
>