Subroutine euler22m(d1x,d2x,d1y,d2y,ti,tf,xi,xf,yi,yf)
c=================================================
c euler22m.f: Solution of the second-order 2D ODE
c method:     modified Euler (predictor-corrector)
c written by: Alex Godunov
c last revision: 21 October 2006
c-------------------------------------------------
c input ...
c d1x(t,x,y)- function dx/dt   (supplied by a user)
c d2x(t,x,y)- function d2x/dt2 (supplied by a user)
c d1y(t,x,y)- function dy/dt   (supplied by a user)
c d2y(t,x,y)- function d2y/dt2 (supplied by a user)
c    where x(2) and y(2) (x(1)-position, x(2)-speed, etc.)
c ti 	- initial time
c tf  - time for a solution
c xi(2) - initial position and speed for x component
c yi(2) - initial position and speed for y component
c
c output ...
c xf(2) - solutions (x position and speed) at point tf
c yf(2) - solutions (y position and speed) at point tf
c=================================================
      implicit none
      Real*8 d1x, d2x, d1y, d2y, ti, tf
      Real*8 xi(2), xf(2), yi(2), yf(2)
      Real*8 h,t, x1, x2, y1, y2
      Real*8 k1x(2),k2x(2),k3x(2),k4x(2),k1y(2),k2y(2),k3y(2),k4y(2)
      h = tf-ti
      t = ti
c*** Euler
      xf(1) = xi(1) + h*d1x(t,xi,yi)
      xf(2) = xi(2) + h*d2x(t,xi,yi)
      yf(1) = yi(1) + h*d1y(t,xi,yi)
      yf(2) = yi(2) + h*d2y(t,xi,yi)
c*** modified Euler
      xf(1) = xi(1) + (d1x(t,xi,yi)+d1x(t,xf,yf))*0.5*h
      xf(2) = xi(2) + (d2x(t,xi,yi)+d2x(t,xf,yf))*0.5*h
      yf(1) = yi(1) + (d1y(t,xi,yi)+d1y(t,xf,yf))*0.5*h
      yf(2) = yi(2) + (d2y(t,xi,yi)+d2y(t,xf,yf))*0.5*h
      Return
      End