Subroutine rk4_2d(d1x,d2x,ti,xi,vi,tf,xf,vf)
c=================================================
c rk4_2d.f:  Solution of the second-order 1D ODE
c method:    Runge-Kutta 4th-order
c written by: Alex Godunov (October 2006)
c-------------------------------------------------
c input ...
c d1x(t,x,v)- function dx/dt   (supplied by a user)
c d2x(t,x,v)- function d2x/dt2 (supplied by a user)
c ti 	- initial time
c xi  - initial position
c vi  - initial velocity (first order derivative)
c tf  - time for a solution
c
c output ...
c xf  - solution at point tf
c vf  - velocity at point tf
c=================================================
      Real*8 d1x,d2x,ti,xi,vi,tf,xf,vf
      Real*8 h,t,k1x,k2x,k3x,k4x,k1v,k2v,k3v,k4v

      h = tf-ti
      t = ti

      k1x = h*d1x(t,xi,vi)
      k1v = h*d2x(t,xi,vi)

      k2x = h*d1x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0)
      k2v = h*d2x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0)

      k3x = h*d1x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0)
      k3v = h*d2x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0)

      k4x = h*d1x(t+h,xi+k3x,vi+k3v)
      k4v = h*d2x(t+h,xi+k3x,vi+k3v)

      xf = xi + (k1x + 2.0*(k2x+k3x) + k4x)/6.0
      vf = vi + (k1v + 2.0*(k2v+k3v) + k4v)/6.0

      Return
      End