3rd Order Runge-Kutta Matlab I'm trying to create a Matlab function to use a matrix form of the 3rd order Runge-Kutta algorithm. I have working code to use the standard RK3 algorithm but I'm struggling to understand how to handle a system of equations. Here is the exact question:

Here is the information I was provided, but I don't understand how to find b or how to implement it as a function handle in Matlab. Can someone give me a few pointers please?


 A: The cited method is slightly wrong. The correct method has the Butcher tableau (see, for instance, the slide of 3rd order methods in https://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/low-order-RK.pdf)
\begin{array}{c|ccc}
0\\
1&1\\
\frac12&\frac14&\frac14\\
\hline
&\frac16&\frac16&\frac23
\end{array}
which can be implemented (quite redundantly) as 
\begin{align}
&&&&k_1&=f(x,y)\\  
y^{(1)}&=y+hk_1&&=y+hf(x,y)& k_2&=f(x+h,y^{(1)})\\
y^{(2)}&=y+\tfrac14hk_1+\tfrac14hk_2&&=\tfrac34y+\tfrac14y^{(1)}+\tfrac14hf(x+h,y^{(1)})& k_3&=f(x+\tfrac12h,y^{(2)})\\
y_+&=y+h(\tfrac16k_1+\tfrac16k_2+\tfrac23k_3)&&=\tfrac13y+\tfrac23y^{(2)}+\tfrac23hf(x+\tfrac12h,y^{(2)})
\end{align}
Thus in the formula for $y_{n+1}$ there is a factor $\frac12$ missing in $b(x_n+\frac12h)$. This error will reduce the order of the method, most likely to order $1$ in the case where the inhomogeneity $b$ is not constant.

As to your question, you define
function b = bvector(x)
    // b components = functions of x
end

(1/20/17, moved up from comments from 11/17/15) In the first order linear ODE $y'(x)−Ay(x)=b(x)$, the vector b(x)
is the inhomogeneity. As the example is homogeneous, you simply get b=[ 0; 0]. It is there just to have more generality in the problem class you can solve.

Final comment: Stability demands that $λh∈[−2.51,0]$
for all (real) eigenvalues $λ$ of $A$, the stability region in the left half plane of the complex plane is more complicated than just a circle or rectangle over that interval, but that should give an idea. In the example, this severely restricts the step size as there is one eigenvalue $λ=-1000$.
