I will outline the process and you can fill in the calculations.
We have our system as:
$$
\left\{\begin{array}{l}
\frac{dy}{dx} = z \\
\frac{dz}{dx} = 6y - z
\end{array}\right.
$$
With $y(0)=3$ and $z(0)=1$.
We must do the calculations in a certain order as there are dependencies between the numerical calculations. This order is:
- $k_0 = hf(x_i,y_i,z_i)$
$l_0 = hg(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)$
$l_1 = hg(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)$
$l_2 = hg(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)$
$l_3 = hg(x_i+h,y_i+k_2,z_i+l_2)$
$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)$
We typically need some inputs for the algorithm:
- A range that we want to do the calculations over: $a \le t \le b$, lets use $a = 0, b = 1$.
- The number of steps $N$, say $N = 10$.
- The steps size $h = \dfrac{b-a}{N} = \dfrac{1}{10}$
The system we are solving is:
$$ \frac{dy}{dx} = f(x,y,z) = z \\ \frac{dz}{dx}=g(x,y,z) = 6y - z$$
Doing the calculations using the above order for the first time step $i= 0, t_0 = 0 = x_0$, yields:
- $k_0 = hf(x_0,y_0,z_0) = \dfrac{1}{10}(z_0) = \dfrac{1}{10}(1) = \dfrac{1}{10}$
$l_0 = hg(x_0,y_0,z_0) = \dfrac{1}{10}(6y_0 - z_0) = \dfrac{1}{10}(6 \times 3 - 1) = \dfrac{1}{10}(17)$
$k_1 = hf(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_0,z_0+\frac{1}{2}l_0) = \dfrac{1}{10}(1 + \dfrac{1}{2}\dfrac{1}{10}(17)) ~~$(You please continue the calcs.)
$l_1 = hg(x_0+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_0+\frac{1}{2}l_0)$
$k_2 = hf(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_1,z_0+\frac{1}{2}l_1)$
$l_2 = hg(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_1,z_0+\frac{1}{2}l_1)$
$k_3 = hf(x_0+h,y_0+k_2,z_0+l_2)$
$l_3 = hg(x_0+h,y_0+k_2,z_0+l_2)$
$y_{1}=y_0 + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$
- $z_{1}=z_0 + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$
You now have $x_1$ and $z_1$ which you need for the next time step after all of the intermediate (in order again). Now, we move on to the next time step $i = 1, t_1 = t_0 + h = \dfrac{1}{10} = x_1$, so we have:
- $k_0 = hf(x_1,y_1,z_1) = \dfrac{1}{10}(z_1)$
$l_0 = hg(x_1,y_1,z_1) = \dfrac{1}{10}(6y_1 - z_1)$
$k_1 = hf(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_0,z_1+\frac{1}{2}l_0)$
$l_1 = hg(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_0,z_1+\frac{1}{2}l_0)$
$k_2 = hf(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_1,z_1+\frac{1}{2}l_1)$
$l_2 = hg(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_1,z_1+\frac{1}{2}l_1)$
$k_3 = hf(x_1+h,y_1+k_2,z_1+l_2)$
$l_3 = hg(x_1+h,y_1+k_2,z_1+l_2)$
$y_{2}=y_1 + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$
- $z_{2}=z_1 + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$
Continue this for $10$ time steps. Your final result should match closely (assuming the numerical algorithm is stable for this problem) to the exact solution. You will compare $z_{10}$ to the exact result. The exact solution is:
$$y(x) = e^{-3 x}+2 e^{2 x}$$
If we find $y(1) = \dfrac{1}{e^3} + 2 e^2 = 14.8278992662291643974401973...$.
Assuming you're using method of lines.
Let the original initial-boundary problem be
$$
u_t = u_{xx}\\
u(0, x) = f(x)\\
u(t, 0) = a(t)\\
u(t, 1) = b(t).
$$
Introduce a set of points $x_j = jh,\; j = 0,1,\dots, N,\;Nh = 1$.
Let $u_j(t) = u(t, x_j)$. Note that $u_0(t) = a(t),\; u_N(t) = b(t)$ are already known. Unknown are $u_j(t),\; j = 1, 2, \dots, N-1$.
Then $u_{xx}(t, x_j)$ can be approximated as
$$
u_{xx}(t, x_j) \approx \frac{u_{j-1}(t) - 2u_j(t) + u_{j+1}(t)}{h^2}.
$$
Plugging that into PDE gives
$$
u'_j(t) = \frac{u_{j-1}(t) - 2u_j(t) + u_{j+1}(t)}{h^2}
$$
a system of $N-1$ ODEs with initial conditions $u_j(0) = f(x_j)$. That could be solved using any RK method (provided that method is stable). For explicit Euler method that would be
$$
\frac{u_j(t_{n+1}) - u_j(t_n)}{\tau} = \frac{u_{j-1}(t_n) - 2u_j(t_n) + u_{j+1}(t_n)}{h^2}, \; j = 1, 2, \dots, N-1\\
u_j(0) = f(t_j)\\
u_0(t) = a(t), \; u_N(t) = b(t).
$$
This well known explicit scheme is stable when $\frac{\tau}{h^2} \leqslant \frac{1}{2}$.
Edit. For those who ask how to use this method with higher order RK methods. Let's take for example an RK2 method (explicit midpoint) with the following Butcher's tableau:
$$
\begin{array}{c|cc}
0 & 0 & 0\\
1/2 & 1/2 & 0\\
\hline
& 0 & 1
\end{array}
$$
Applied to ODE in form $\mathbf u' = \mathbf F(t, \mathbf u)$ this method expands to
$$
\mathbf r = \mathbf F(t_n, \mathbf u_n)\\
\mathbf s = \mathbf F\left(t_n + \frac{\tau}{2}, \mathbf u_n + \frac{\tau}{2} \mathbf r\right)\\
\frac{\mathbf u_{n+1} - \mathbf u_n}{\tau} = \mathbf s
$$
I've used $\mathbf r$ and $\mathbf s$ instead of common $\mathbf k_{1,2}$ to reduce the number of indices involved. Here $\mathbf r$ and $\mathbf s$ are intermediate values that depend solely on values of $u_j$ at $t = t_n$.
For our case the right hand side of the ODE is
$$
F_j(t, u_1, u_2, \dots, u_{N-1}) = \frac{1}{h^2}\begin{cases}
u_0(t) - 2 u_1 + u_2, &j = 1\\
u_{j-1} - 2 u_j + u_{j+1}, &1 < j < N-1\\
u_{N-2} - 2 u_{N-1} + u_N(t), &j = N-1
\end{cases}
$$
Note that $u_0(t)$ and $u_N(t)$ are given explicitly as $a(t)$ and $b(t)$. This is why I have separated cases $j=1$ and $j = N-1$ in the definition of $F_j$.
Putting this altogether gives
$$
r_j = \frac{u_{j-1}(t_n) - 2 u_j(t_n) + u_{j-1}(t_n)}{h^2}\\
s_j = \frac{1}{h^2}\begin{cases}
u_0\left(t + \frac{\tau}{2}\right) - 2 \left(u_{1}(t_n) + \frac{\tau}{2}r_{1}\right) + \left(u_2(t_n) + \frac{\tau}{2}r_2\right), &j = 1\\
\left(u_{j-1}(t_n) + \frac{\tau}{2}r_{j-1}\right) - 2 \left(u_j(t_n) + \frac{\tau}{2}r_j\right) + \left(u_{j+1}(t_n) + \frac{\tau}{2}r_{j+1}\right), &1 < j < N-1\\
\left(u_{N-2}(t_n) + \frac{\tau}{2}r_{N-2}\right) - 2 \left(u_{N-1}(t_n) + \frac{\tau}{2}r_{N-1}\right) + u_N\left(t + \frac{\tau}{2}\right), &j = N-1
\end{cases}\\
\frac{u_j(t_{n+1}) - u_j(t_n)}{\tau} = s_j
$$
Note that $r_j$ and $s_j$ are helper values to step from $u_j(t_n)$ to $u_j(t_{n+1})$ and are different for each time step. One may want to attribute values $r_j$ and $s_j$ to some moment of time. A reasonable choice would be attributing all each value $\mathbf k_i$ with moment $t_n + \tau c_i$. Here $c_i$ is the first column of the Butcher's tableau.
$$
r_j(t_n) = \frac{u_{j-1}(t_n) - 2 u_j(t_n) + u_{j-1}(t_n)}{h^2}\\
s_j\left(t_n + \frac{\tau}{2}\right) = \frac{1}{h^2} \times \\ \times \begin{cases}
u_0\left(t + \frac{\tau}{2}\right) - 2 \left(u_{1}(t_n) + \frac{\tau}{2}r_{1}(t_n)\right) + \left(u_2(t_n) + \frac{\tau}{2}r_2(t_n)\right), &j = 1\\
\left(u_{j-1}(t_n) + \frac{\tau}{2}r_{j-1}(t_n)\right) - 2 \left(u_j(t_n) + \frac{\tau}{2}r_j(t_n)\right) + \left(u_{j+1}(t_n) + \frac{\tau}{2}r_{j+1}(t_n)\right), &1 < j < N-1\\
\left(u_{N-2}(t_n) + \frac{\tau}{2}r_{N-2}(t_n)\right) - 2 \left(u_{N-1}(t_n) + \frac{\tau}{2}r_{N-1}(t_n)\right) + u_N\left(t + \frac{\tau}{2}\right), &j = N-1
\end{cases}\\
\frac{u_j(t_{n+1}) - u_j(t_n)}{\tau} = s_j\left(t_n + \frac{\tau}{2}\right)
$$
While this is the answer to the question "How to apply RK2 to this ODE" I really don't like the final form. Instead let's write the same method in a slightly different form:
$$
\frac{\mathbf u_{n+1/2} - \mathbf u_n}{\tau / 2} = \mathbf F(t_n, \mathbf u_n)\\
\frac{\mathbf u_{n+1} - \mathbf u_n}{\tau} = \mathbf F(t_n + \frac{\tau}{2}, \mathbf u_{n+1/2}).
$$
One can check that the methods are the same by plugging $\mathbf u_{n+1/2} = \mathbf u_n + \frac{\tau}{2} \mathbf r$. Just like $\mathbf r$ and $\mathbf s$ the $\mathbf u_{n+1/2}$ is a helper value used to perform a step from $\mathbf u_n$ to $\mathbf u_{n+1}$.
Applied to our ODE this method gives
$$
\frac{u_j(t_{n+1/2}) - u_j(t_n)}{\tau / 2} = \frac{u_{j-1}(t_n) - 2 u_j(t_n) + u_{j-1}(t_n)}{h^2}, \quad j = 1, 2, \dots, N-1\\
u_0(t_{n+1/2}) = a\left(t_n + \frac{\tau}{2}\right), \quad
u_N(t_{n+1/2}) = b\left(t_n + \frac{\tau}{2}\right),\\
\frac{u_j(t_{n+1}) - u_j(t_n)}{\tau / 2} = \frac{u_{j-1}(t_{n+1/2}) - 2 u_j(t_{n+1/2}) + u_{j-1}(t_{n+1/2})}{h^2}, \quad j = 1, 2, \dots, N-1\\
u_0(t_{n+1}) = a\left(t_n + \tau\right), \quad
u_N(t_{n+1}) = b\left(t_n + \tau\right).
$$
Though not strictly necessary I have defined values $u_0(t_{n+1/2})$ and $u_N(t_{n+1/2})$ to get rid of treating $j=1$ and $j=N-1$ as separate cases.
Best Answer
If you have procedures
A(t,x)andb(t,x)that return the system matrix and rhs vector (this is just as placeholders, one can construct these also in-place), then you can do (using python syntax)and use the usual RK4 step
etc.
Similarly in matlab
and the usual RK4 step
etc. There is also an
ode4procedure available in the Matlab studios (? Mathworks side with examples and tutorial videos), but it is for scalar examples, needs to be adapted for systems, mainly in growing the results array.