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.

A Fortran code shown below:

produces the following result

!Runge-Kutta Fourth Order Method

!For 2nd Order Differentiation Equation

!First you have to define the function

F(x,y,z) = z !dy/dx

G(x,y,z) = 6*y-z !dz/dx = d2y/dx2

INTEGER :: n,i

REAL :: k1,l1,k2,l2,k3,l3,k4,l4    !Most Important

Write (*,*) "Given Equation '(y2)-6(y1)+(y0)=0'"

Write (*,*) "Xo=0, Yo=3, Zo=Y'o=1, Xn=1, n=?"

Xo=0    !Given Condition

Yo=3    !Given Condition

Zo=1    !Given Condition

Xn=1    !Given Condition

read (*,*) n    !n=number of Intercept

h=(Xn-Xo)/n

do i=1,n    !you have to do the Calculation 'n' times

k1 = h*F(Xo,Yo,Zo)

l1 = h*G(Xo,Yo,Zo)

k2 = h*F(Xo+h/2,Yo+k1/2,Zo+l1/2)

l2 = h*G(Xo+h/2,Yo+k1/2,Zo+l1/2)

k3 = h*F(Xo+h/2,Yo+k2/2,Zo+l2/2)

l3 = h*G(Xo+h/2,Yo+k2/2,Zo+l2/2)

k4 = h*F(Xo+h,Yo+k3,Zo+l3)

l4 = h*G(Xo+h,Yo+k3,Zo+l3)

!Sum Up

Yn = Yo+(k1+2*k2+2*k3+k4)/6

Zn = Zo+(l1+2*l2+2*l3+l4)/6

!Operation for Next calculation

Xo=Xo+h     !(+h) than previous Term

Yo=Yn    !Now Yn becomes Yo

Zo=Zn    !Now Zn becomes Zo

End Do

Write (*,*) "Xn,Yn =",Xo,Yo

Stop

End

