!3rd order Runge-Kutta algorithm routine for a system of 
!differential equations.

subroutine RK_general(f, x0, y0, h, x1, y1, d)
implicit none

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

do i = 1, d
   k0(i)=f(x0, y0(1), y0(2), i)
enddo

call RK_k_half_general(f, x0, y0, h, k0, k_half, d)
call RK_k_one_general(f, x0, y0, k0, h, k1, d)

do i = 1, d
   k(i) = (1.0/6.0)*k0(i) + (4.0/6.0)*k_half(i) + (1.0/6.0)*k1(i)
   y1(i) = y0(i) + k(i)*h
enddo

x1 = x0 + h
end subroutine RK_general