%% Numerical solution of initial value problem using ode45 and Euler method
%% Using ode45

clearvars; clf

% We want to find the solution of the IVP
%        y'=y-t,\quad y(0)=\frac12
% for 0\le t \le T with T=1

f = @(t,y) y-t
t0 = 0; y0 = .5; T = 1;
dirfield(f,0:.1:1,.5:.1:1); hold on

% We find the "exact" solution using ode45:

[ts,ys] = ode45(f,[t0,T],y0);
yex = ys(end)
plot(ts,ys,'b')
text(T,yex,'        exact')
xlabel('t'); title('Numerical solution of IVP using ode45')

%
%% Euler method
% The function EulerMethod is defined at the end of this file.

for n = 2.^(0:6)
  [ts,ys] = EulerMethod(f,[t0,T],y0,n);
  plot(ts,ys,'.-') 
  Y = ys(end);          % approximation for y(T)
  fprintf('n = %2g, y(T)≈%8g, error = %8g \n',n,Y,Y-yex)
  if n<=4
    disp([ts,ys])
  end
  text(T,Y,sprintf('        n=%g',n))
end
hold off
title('Euler Method using 1,2,4,...,64 steps')

%
%% Code for Euler method

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;
ts = zeros(n+1,1); ys = zeros(n+1,1);   % vectors for t-values, y-values
ts(1) = t; ys(1) = y;
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 vectors ts, ys
end
end