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

subroutine RK_k_one_general(f, x0, y0, k0, h, k1, d)
implicit none

     integer :: d, i
     real :: y0(d), k0(d), k1(d), x0, x1, h
     real, external :: f

x1 = x0 + h 

do i = 1, d
   k1(i) = f(x1, y0(1)+h*k0(i), y0(2)+h*k0(i), i)
enddo

end subroutine RK_k_one_general