function dx=f(x,t)
  r=sqrt(x(1)**2+x(2)**2);
  dx(1)=x(3);
  dx(2)=x(4);
  dx(3)=-x(1)/r**3;
  dx(4)=-x(2)/r**3;