If the IVP is linear, that is, the function **f** may be written as

**f**(t,

**u**) =

**Mu**+

**b**

for a matrix **M** and vector **b**, then we can very easily
use Matlab to solve such systems.

For example, given the IVP

^{(1)}(t) = -0.5 x(t) + 1.9 y(t) + 0.4

y

^{(1)}(t) = -2.6 x(t) − 0.7 y(t) + 0.8

x(0) = 1

y(0) = 1

then we may assign **M** and **b**:

>> M = [-0.5 1.9; -2.6 0.7]; >> b = [0.4 0.8]';

We will store the iterations in the 2 × *N* matrix and
assign the initial values to be the first column:

>> N = 10000; >> x = zeros( 2, N ); >> x(:, 1) = [1 1]';

Next, we iterate Euler's method by assigning to the next column:

>> h = 0.01; >> for i=2:N x(:, i) = x(:, i - 1) + h*(M*x(:, i - 1) + b); end;

To plot this data, we pass `plot` the first row
(the *x* values) and the second row (the *y* values):

>> plot( x(1, :), x(2, :) )

The result of the first 10000 iterations is shown in Figure 1.

Figure 1. The first 10000 approximations of the solution to the system of IVPs (click to enlarge).

You may note that the solution appears to converge to a point.
That point is the value of **u** such that **Mu** + **b** = **0**,
or **Mu** = −**b**. Thus, to find this point, we solve:

>> M \ (-b) ans = 0.340264650283554 -0.120982986767486

# Heun's Method

To implement Heun's method, we need only one change:

>> for i=2:N K0 = M*x(:, i - 1) + b; K1 = M*(x(:, i - 1) + h*K0) + b; x(:, i) = x(:, i - 1) + h*(K0 + K1)/2; end;

# 4th-order Runge Kutta Method

To implement the 4th-order Runge Kutta method, we need only one change:

>> for i=2:N K0 = M*x(:, i - 1) + b; K1 = M*(x(:, i - 1) + h/2*K0) + b; K2 = M*(x(:, i - 1) + h/2*K1) + b; K3 = M*(x(:, i - 1) + h*K2) + b; x(:, i) = x(:, i - 1) + h*(K0 + 2*K1 + 2*K2 + K3)/6; end;

Copyright ©2006 by Douglas Wilhelm Harder. All rights reserved.