Although this answer contains the same content as Amzoti's answer, I think it's worthwhile to see it another way.
In general consider if you had $m$ first-order ODE's (after appropriate decomposition). The system looks like
\begin{align*}
\frac{d y_1}{d x} &= f_1(x, y_1, \ldots, y_m) \\
\frac{d y_2}{d x} &= f_2(x, y_1, \ldots, y_m) \\
&\,\,\,\vdots\\
\frac{d y_m}{d x} &= f_m(x, y_1, \ldots, y_m) \\
\end{align*}
Define the vectors $\vec{Y} = (y_1, \ldots, y_m)$ and $\vec{f} = (f_1, \ldots, f_m)$, then we can write the system as
$$\frac{d}{dx} \vec{Y} = \vec{f}(x,\vec{Y})$$
Now we can generalize the RK method by defining
\begin{align*}
\vec{k}_1 &= h\vec{f}\left(x_n,\vec{Y}(x_n)\right)\\
\vec{k}_2 &= h\vec{f}\left(x_n + \tfrac{1}{2}h,\vec{Y}(x_n) + \tfrac{1}{2}\vec{k}_1\right)\\
\vec{k}_3 &= h\vec{f}\left(x_n + \tfrac{1}{2}h,\vec{Y}(x_n) + \tfrac{1}{2}\vec{k}_2\right)\\
\vec{k}_4 &= h\vec{f}\left(x_n + h, \vec{Y}(x_n) + \vec{k}_3\right)
\end{align*}
and the solutions are then given by
$$\vec{Y}(x_{n+1}) = \vec{Y}(x_n) + \tfrac{1}{6}\left(\vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4\right)$$
with $m$ initial conditions specified by $\vec{Y}(x_0)$. When writing a code to implement this one can simply use arrays, and write a function to compute $\vec{f}(x,\vec{Y})$
For the example provided, we have $\vec{Y} = (y,z)$ and $\vec{f} = (z, 6y-z)$. Here's an example in Fortran90:
program RK4
implicit none
integer , parameter :: dp = kind(0.d0)
integer , parameter :: m = 2 ! order of ODE
real(dp) :: Y(m)
real(dp) :: a, b, x, h
integer :: N, i
! Number of steps
N = 10
! initial x
a = 0
x = a
! final x
b = 1
! step size
h = (b-a)/N
! initial conditions
Y(1) = 3 ! y(0)
Y(2) = 1 ! y'(0)
! iterate N times
do i = 1,N
Y = iterate(x, Y)
x = x + h
end do
print*, Y
contains
! function f computes the vector f
function f(x, Yvec) result (fvec)
real(dp) :: x
real(dp) :: Yvec(m), fvec(m)
fvec(1) = Yvec(2) !z
fvec(2) = 6*Yvec(1) - Yvec(2) !6y-z
end function
! function iterate computes Y(t_n+1)
function iterate(x, Y_n) result (Y_nplus1)
real(dp) :: x
real(dp) :: Y_n(m), Y_nplus1(m)
real(dp) :: k1(m), k2(m), k3(m), k4(m)
k1 = h*f(x, Y_n)
k2 = h*f(x + h/2, Y_n + k1/2)
k3 = h*f(x + h/2, Y_n + k2/2)
k4 = h*f(x + h, Y_n + k3)
Y_nplus1 = Y_n + (k1 + 2*k2 + 2*k3 + k4)/6
end function
end program
This can be applied to any set of $m$ first order ODE's, just change m in the code and change the function f to whatever is appropriate for the system of interest. Running this code as-is yields
$ 14.827578509968953 \qquad 29.406156886687729$
The first value is $y(1)$, the second $z(1)$, correct to the third decimal point with only ten steps.
