%% Numerical solution of IVP for 2nd order ODE

clearvars; clf

% We want to find the solution of the IVP
%        y''+3y'+2y=t,
% \qquad y(0)=2, \quad y'(0)=-1
% for t\in [0,5].
% This is a second order ODE.
% We first have to rewrite this as a first order system:
% Let y_1(t):=y(t),\quad y_2(t):=y'(t) with the vector valued function \vec y(t) = \left[
% \matrix{y_1(t) \cr y_2(t)}
% \right].  Note that
%         y_1'=y' = y_2,  \qquad y_2' = y'' = t-3y'-2y = t-3y_2-2y_1
% We obtain the first order system
%            \vec y' = \vec f(t,\vec y),\qquad \vec y(t_0)=\vec y^{(0)}
% where we have
%            \vec y' = \left[
% \matrix{y_2 \cr t-3y_2-2y_1}
% \right],
% \qquad \vec y(t_0)=\left[
% \matrix{2 \cr -1}
% \right]
% We define the function \vec f(t,\vec y) 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=\frac12

h = 1/2;

% We evaluate \vec s = \vec f(t_0,\vec y^{(0)}):

s = f(0,y0)

% and obtain for t_1=\frac12 the approximation \vec y^{(1)} :=\vec y^{(0)} + h\cdot\vec s

y1 = y0 + h*s

% We obtain the approximations y(\frac12) \approx 1.5, \quad
% y'(\frac12) \approx -1.5.
%% Use 10 steps of the Euler method for t\in[0,5]
% Here the step size is h = \frac{T-t_0}n = \frac12.
% The function EulerMethod is defined at the end of this file.

[ts,ys] = EulerMethod(f,[0,5],y0,10);
disp([ts,ys])

% Note that ys has two columns containing the values of y(t) and y'(t).
% Plot the Euler approximation for y(t):

plot(ts,ys(:,1),'b.-'); 
xlabel('t'); title('Euler approximation for IVP');hold on

%
%% Using ode45
% Find the solution \vec y(t) for t\in [0,5]:

[ts,ys] = ode45(f,[0,5],y0);
disp([ts,ys])   % shows t, y(t), y'(t) in each row (scroll to see all rows)

% Note that ys has two columns containing the values of y(t) and y'(t).
% Show the value for y(T):

yex = ys(end,1)

% Plot the solution y(t):

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 \vec y^ 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