!Runge-Kutta method for integration of ODE's.
!Subroutine returning k1 = f(x1,y0+h*k0).

subroutine RK_k_one(f, x0, y0, h, k1, x1)
implicit none

     real :: y0, k0, k1, x0, x1, h
     real, external :: f

k0=f(x0,y0)
x1=x0 + h ; k1=f(x1,y0+h*k0)

end subroutine RK_k_one