Numerical solution of IVP for 2nd order ODE

clearvars; clf
We want to find the solution of the IVP
for .
This is a second order ODE.
We first have to rewrite this as a first order system:
Let with the vector valued function . Note that
We obtain the first order system
where we have
We define the function in Matlab:
f = @(t,y) [ y(2) ; t-3*y(2)-2*y(1) ];
Our initial vector is
y0 = [2;-1];

Perform one step of Euler method with

h = 1/2;
We evaluate :
s = f(0,y0)
s = 2×1
-1 -1
and obtain for the approximation
y1 = y0 + h*s
y1 = 2×1
1.5000 -1.5000
We obtain the approximations .

Use 10 steps of the Euler method for

Here the step size is .
The function EulerMethod is defined at the end of this file.
[ts,ys] = EulerMethod(f,[0,5],y0,10);
disp([ts,ys])
0 2.0000 -1.0000 0.5000 1.5000 -1.5000 1.0000 0.7500 -0.5000 1.5000 0.5000 0 2.0000 0.5000 0.2500 2.5000 0.6250 0.3750 3.0000 0.8125 0.4375 3.5000 1.0312 0.4688 4.0000 1.2656 0.4844 4.5000 1.5078 0.4922 5.0000 1.7539 0.4961
Note that ys has two columns containing the values of and .
Plot the Euler approximation for :
plot(ts,ys(:,1),'b.-');
xlabel('t'); title('Euler approximation for IVP');hold on

Using ode45

Find the solution for :
[ts,ys] = ode45(f,[0,5],y0);
disp([ts,ys]) % shows t, y(t), y'(t) in each row (scroll to see all rows)
0 2.0000 -1.0000 0.0502 1.9486 -1.0430 0.1005 1.8954 -1.0728 0.1507 1.8410 -1.0910 0.2010 1.7860 -1.0992 0.2876 1.6909 -1.0939 0.3742 1.5972 -1.0688 0.4608 1.5062 -1.0285 0.5474 1.4192 -0.9772 0.6394 1.3322 -0.9145 0.7314 1.2512 -0.8461 0.8234 1.1766 -0.7741 0.9154 1.1087 -0.7006 1.0188 1.0406 -0.6183 1.1222 0.9809 -0.5374 1.2256 0.9294 -0.4589 1.3289 0.8858 -0.3837 1.4410 0.8472 -0.3067 1.5531 0.8169 -0.2344 1.6653 0.7945 -0.1671 1.7774 0.7792 -0.1048 1.8955 0.7705 -0.0445 2.0136 0.7685 0.0105 2.1317 0.7728 0.0607 2.2499 0.7827 0.1062 2.3612 0.7967 0.1450 2.4726 0.8149 0.1804 2.5840 0.8368 0.2124 2.6954 0.8620 0.2413 2.8204 0.8941 0.2706 2.9454 0.9296 0.2966 3.0704 0.9681 0.3198 3.1954 1.0094 0.3404 3.3204 1.0531 0.3587 3.4454 1.0990 0.3750 3.5704 1.1468 0.3894 3.6954 1.1963 0.4022 3.8204 1.2473 0.4135 3.9454 1.2996 0.4236 4.0704 1.3531 0.4324 4.1954 1.4077 0.4403 4.3204 1.4631 0.4473 4.4454 1.5194 0.4534 4.5704 1.5765 0.4589 4.6954 1.6341 0.4637 4.7715 1.6695 0.4663 4.8477 1.7052 0.4688 4.9238 1.7409 0.4710 5.0000 1.7769 0.4732
Note that ys has two columns containing the values of and .
Show the value for :
yex = ys(end,1)
yex = 1.7769
Plot the solution :
plot(ts,ys(:,1),'k'); hold off
title('approximations from ode45 (black) and Euler method (blue)')

Code for Euler method

We modify the code from the scalar case so that it works for systems: for each time we store the vector as a row of the array ys.
function [ts,ys] = EulerMethod(f,interval,y,n)
% [ts,ys] = EulerMethod(f,[t0,T],y0,n)
% initial value problem y'=f(t,y), y(t0) = y0
% use n steps of Euler method to approximate solution in interval [t0,T]
t = interval(1); T = interval(2);
h = (T-t)/n;
m = size(y,1); % size of vector y
ts = zeros(n+1,1); ys = zeros(n+1,m); % initialize vector for t-values, array for y-values
ts(1) = t; ys(1,:) = y.'; % y.' is transpose, gives row vector
for k=1:n
s = f(t,y);
t = t + h;
y = y + s*h;
ts(k+1) = t; ys(k+1,:) = y.'; % store t,y in vector ts, array ys
end
end