subroutine orbit (t, y, yp)
      double precision t, y(4), yp(4), r, alfasq
      double precision dsqrt
      common alfasq
      r = y(1)*y(1) + y(2)*y(2)
      r = r*dsqrt(r)/alfasq
      yp(1) = y(3)
      yp(2) = y(4)
      yp(3) = -y(1)/r
      yp(4) = -y(2)/r
      return
      end

.