Help with using the Runge-Kutta 4th order method on a system of 2 first order ODE's. The original ODE I had was $$ \frac{d^2y}{dx^2}+\frac{dy}{dx}-6y=0$$ with $y(0)=3$ and $y'(0)=1$. Now I can solve this by hand and obtain that $y(1) = 14.82789927$. However I wish to use the 4th order Runge-Kutta method, so I have the system:
$$
\left\{\begin{array}{l}
  \frac{dy}{dx} = z \\
  \frac{dz}{dx} = 6y - z
\end{array}\right.
$$
With $y(0)=3$ and $z(0)=1$.
Now I know that for two general 1st order ODE's $$ \frac{dy}{dx} = f(x,y,z) \\ \frac{dz}{dx}=g(x,y,z)$$ The 4th order Runge-Kutta formula's for a system of 2 ODE's are: $$ y_{i+1}=y_i + \frac{1}{6}(k_0+2k_1+2k_2+k_3) \\ z_{i+1}=z_i + \frac{1}{6}(l_0+2l_1+2l_2+l_3) $$ Where $$k_0 = hf(x_i,y_i,z_i) \\ k_1 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0) \\ k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1) \\ k_3 = hf(x_i+h,y_i+k_2,z_i+l_2) $$ and $$l_0 = hg(x_i,y_i,z_i) \\ l_1 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0) \\ l_2 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1) \\ l_3 = hg(x_i+h,y_i+k_2,z_i+l_2)$$
My problem is I am struggling to apply this method to my system of ODE's so that I can program a method that can solve any system of 2 first order ODE's using the formulas above, I would like for someone to please run through one step of the method, so I can understand it better.

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.
