Numerical solution of an ODE system of equations using RK4 I have given an assignment to find the solution to the ODE system of equations as follow:
$$\begin{cases} x_1' = x_1 + x_2 \\ x_2' = -3x_1 -10x_2 + x_2 ^2\end{cases}$$
With initial conditions: $x_1(0) = 1, \quad x_2(0)=10$ on interval $[0,10]$
The above system of questions differs in two aspects from the normal questions of my textbook so I am confused and don't know how to program the algorithm in Matlab:


*

*It has two intermingled equations. In contrary to textbook samples which has only one in form of: $x'=f(t,x)$

*$x_1'$ and $x_2'$ don't have any $t$ variable! which makes me even more confused.


I appreciate it if you elaborate  on the above two points. Thanks.
------------------------------------------Part2-----------------------------------------
Ok with above clarifications I started to program in Matlab, after running my code, $x_1(t)$ and $x_2(t)$ seems to change (probably becomes unstable) after changing the step number which changes h respectively.
Here is the code:
function [sol_x1, sol_x2] = rk4 (m)
    a = 0;      % interval = [a, b] = [0, 10]
    b = 10;
    h = (b-a) / m; % m is  number of steps

    Sx1 = zeros (1,m+1);      % Solutions will be saved here
    Sx2 = zeros (1,m+1);

    Sx1(1) = 1;    % initial values
    Sx2(1) = 10;

    for j = 1:m
        x1 = Sx1(j);
        x2 = Sx2(j);

        fprintf('x1 = %i, x2 = %i \n', x1, x2);

        x1_k1 = h * f1(x1, x2);
        x2_k1 = h * f2(x1, x2);

        x1_k2 = h * f1(x1 + x1_k1/2, x2 + x2_k1/2);
        x2_k2 = h * f2(x1 + x1_k1/2, x2 + x2_k1/2);

        x1_k3 = h * f1(x1 + x1_k2/2, x2 + x2_k2/2);
        x2_k3 = h * f2(x1 + x1_k2/2, x2 + x2_k2/2);

        x1_k4 = h * f1(x1 + x1_k3, x2 + x2_k3);
        x2_k4 = h * f2(x1 + x1_k3, x2 + x2_k3);


        next_x1 = x1 + (x1_k1 + (2 * x1_k2) + (2 * x1_k3) + x1_k4) / 6;
        next_x2 = x2 + (x2_k1 + (2 * x2_k2) + (2 * x2_k3) + x2_k4) / 6;

        % no need for t = t + h because x1 and x2 do not depend on t

        %Save the solutions
        fprintf('%i\n', h);
        %display(h)
        Sx1(j+1) = next_x1;
        Sx2(j+1) = next_x2;
    end    

    sol_x1 = Sx1;
    sol_x2 = Sx2;
end

function [y] = f1 (x1, x2) 
    y = x1 + x2;
end

function [y] = f2 (x1, x2) 
    y = (-3 * x1) + (-10 * x2) + (x2 ^ 2);
end

Am I on right path?
$x1(t)$ and $x2(t)$ shoots up to very large numbers pretty quickly. I think something is wrong with my algorithm implementation.
 A: You can think of this system of differential equations as ${\bf x}' = f({\bf x},t)$ where ${\bf x}$ is the vector $\pmatrix{x_1\cr x_2\cr}$, and 
$$ f\left(\pmatrix{x_1\cr x_2\cr},t\right) = \pmatrix{x_1 + x_2\cr −3 x_1−10 x_2+x_2^2}$$
The fact that this doesn't happen to depend on $t$ makes no difference.  The fact that $\bf x$  is a vector rather than a scalar also doesn't matter: the RK4 formulas are exactly the same as in the scalar case.
A: Additionally, you have a problem with the line
h = h + step

while you use h as step size. Since originally step = h you get a rapidly growing step size where it should be constant. What you probably meant was
t = t + h

there is no need for the extra step variable. Also, consider using h=(b-a)/m for greater flexibility, using a single point of constant definition.

Everything else is looking good. You can always test for reliability of the code by varying the step size. The numerical results should stay stable with small variations.
A: Note that the notation $f(x,t)$ is for completeness, in the sense that it is also meant to cover cases where there is no explicit $t$ dependence.
As for the other part of your question, you can write
\begin{align}
x_1' & = f_1(x_1,x_2,t) = \cdots \\
x_2' & = f_2(x_1,x_2,t) = \cdots
\end{align}
If you want me to see if I can find some examples online, or write some pseudo-code, I can do that for you.
