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

subroutine RKSTEP(f, x0, y0, h, d, y1, err)
use MY_STUFF
implicit none

     integer :: i
     integer, intent(in) :: d 
     real(KIND=dbl) :: k0(d), k_half(d), k1(d), k(d), y_midpoint(d), & 
                       & y1_test(d), diff(d)
     real(KIND=dbl), intent(in) :: x0, y0(d), h
     real(KIND=dbl), intent(out) :: y1(d), err(d)
     real, external :: f

!Full step:
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

!2 * half step:
call RK_k_half_general(f, x0, y0, h/2.0, k0, k_half, d)
call RK_k_one_general(f, x0, y0, k0, h/2.0, 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)
   y_midpoint(i) = y0(i) + k(i)*(h/2.0)
enddo

call RK_k_half_general(f, x0+h, y_midpoint, h/2.0, k0, k_half, d)
call RK_k_one_general(f, x0+h, y_midpoint, k0, h/2.0, 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_test(i) = y_midpoint(i) + k(i)*(h/2.0)
enddo

do i = 1, d
   diff(i) = y1_test(i) - y1(i)
   err(i) = ABS(diff(i))/(2.0**3.0 - 1)
enddo

end subroutine RKSTEP