Skip to main content

Section 19.5 Solving systems of differential equations

Conceptually, any of the above methods can be directly applied to systems of \(m\) differential equations, just by interpreting \(y\) as a vector function \(\mathbf y\text{.}\) The notation becomes awkward: now \(y_2\) refer to the second component of \(\mathbf y\text{,}\) while the value \(\mathbf y(t_2)\) might be denoted \(\mathbf y_2\text{.}\) There are also a few changes in Matlab notation regarding y: it will now be a matrix, where rows correspond to components and columns correspond to points of evaluation. The right-hand side function f should return a column vector of the same dimesion \(m\text{.}\)

Solve the system

\begin{align*} y_1' \amp = t - 3y_2\\ y_2' \amp = 2y_1 \end{align*}

with initial condition \(y(0) = \begin{pmatrix}1 \\ 0\end{pmatrix}\text{,}\) on the interval \([0, 10]\) using the trapezoidal method with the step \(h=0.01\text{.}\) Plot the solution.

Answer
f = @(t, y) [t-3*y(2); 2*y(1)]; 
h = 0.01; 
t = 0:h:10;
y0 = [1; 0];
y = y0*ones(size(t));
 
for k = 1:numel(t)-1
    pred = y(:, k) + h*f(t(k), y(:, k));
    y(:, k+1) = y(:, k) + h/2*(f(t(k), y(:, k)) + f(t(k+1), pred));
end

plot(t, y)

The most visible difference is the notation y(:, k) which takes the entire \(k\)th column of the matrix y. Its mathematical meaning is the approximation to \(\mathbf y(t_k)\text{.}\) The initial condition y0 is a column vector, which makes the computation y = y0*ones(size(t)); an outer product: it creates a matrix in which every column is equal to y0. This is convenient because we need the first column to be equal, and the other columns will be rewritten anyway.

The plot displayed by Example 19.5.1 is a time series plot in which one axis represents the independent variable (usually time) and the other axis shows the components of \(\mathbf y\text{.}\) The other option we have is to make a phase plot in which each component of \(\mathbf y\) gets its own axis, and there is no time axis. This is achieved with

plot(y(1,:), y(2,:))

The phase plot may be preferable when the system does not involve time directly (an autonomous system). Try both kinds of plots for the autonomous system

  f = @(t, y) [y(2)-y(2)^3; 2*y(1)];