!Runge-Kutta method for integration of ODE's.
!Subroutine returning k_half = f(x_half,y_half).

subroutine RK_k_half(f, x0, y0, h, k_half)
implicit none

     real :: y0, y_half, k0, k_half, x0, x_half, h
     real, external :: f


k0=f(x0,y0)
y_half=y0 + (h/2.0)*k0 ; x_half=x0 + h/(2.0) ; k_half=f(x_half,y_half)

end subroutine RK_k_half