# 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:

1. It has two intermingled equations. In contrary to textbook samples which has only one in form of: $x'=f(t,x)$
2. $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.

• Presumaby there is a $t$ variable, and $x_1'$ means $\dfrac{dx_1(t)}{dt}$. $\qquad$ Apr 27, 2016 at 17:45
• I agree; I think the OP is confused because the equations are separable (I think that is the term) Apr 27, 2016 at 17:46
• This is true, but that isn't necessary in general. Think of a constant function $f(t)=2$ -- there's no explicit t-dependence, but it is still written "as a function of t". Apr 27, 2016 at 17:48
• Rather than working componentwise as you did here, you can work in vector form, using a function that takes a vector and returns a vector. In this case it would be f(x)=[x(1)+x(2);-3*x(1)-10*x(2)+x(2)^2] You could then write a Matlab implementation exactly like the case for a single ODE because Matlab treats vectors almost exactly the same as it treats numbers.
– Ian
Apr 27, 2016 at 19:54
• Also, I just did the same problem with ode45 and you do indeed get fairly large numbers, on the order of $10^4$ for $x_1$ and $10^2$ for $x_2$
– Ian
Apr 27, 2016 at 19:59

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.

• Your answer cleared my confusion, but I still can't get the code working properly, I deeply appreciate it if you look at my code and see where I have done wrong. Apr 27, 2016 at 18:34

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.

• I think my code is fine. If $m=100$ then $h=0.1$. Therefore, in the for loop $h$ will be increase by '$step = h = 0.1$' variable per each iteration which is correct. I double check by printing $h$ in each iteration and it does not grow rapidly. The problem is somewhere else. After following your correction and adding $x$ it even grows faster! Apr 27, 2016 at 19:17
• Though h=(b-a)/m point is valid, I will add it to the code. Apr 27, 2016 at 19:25
• No, h should be constant, there should be no increase during the loop. Or replace h by step in the computation of the x_k. Apr 27, 2016 at 19:42
• The $+x_2^2$ in the second equation is what is called "positive feedback" and can lead to exploding solutions. The larger $x_2$ gets, the faster it grows. The simple example for this behavior is $u'=u^2$. Apr 27, 2016 at 20:05
• I don't think the rapid growth here is due to the nonlinearity, I think it is mostly just caused by the exponential growth of $x_1$.
– Ian
Apr 28, 2016 at 0:27

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.

• Let me try myself first. I will post it here and then you can fix my errors and then the future users can have a guide to follow when then hit such systems of equations. Apr 27, 2016 at 17:53
• I posted the code, I appreciate if you check. Apr 27, 2016 at 18:33