u1p = u(2);
   u2p = -u(1) / r ^ 3;

   u3p = u(4);
   u4p = -u(3) / r ^ 3;