% BRIEF RECAP of November 5th, 1998 NEEP271 lecture
% written by Fabian Waleffe 98/11/5
%
% ORDINARY DIFFERENTIAL EQUATIONS  (ODEs)
%
% INITIAL VALUE PROBLEMS (IVPs): dy/dt =f(t,y) with y(t_o) = y_o 
%
% Euler's method: March forward by steps of size h
%
%    y(n+1) = y(n) + h f( t_n, y(n) ),  y(0) = y_o
%
% Test case: dy/dt =a y, y(0) = y0 => y(t) = y0 exp(a t)
%

a=2; y0=1;  t0=0; tmax=1;   %define parameters of the problem

% EXACT SOLUTION 

  t=t0+[0:0.02:1]*(tmax-t0); y=y0*exp(a*(t-t0))
  plot(t,y,'m'), xlabel('t'), ylabel('y')
  hold on

% EULER'S METHOD

tns=[t0]; yns=[y0];     % we'll store data in tns and yns

N=4; h=tmax/N; yn=y0;

for i=1:N, 
    yn= yn + h* (a*yn);                 % Euler's method
    tns=[tns; t0+i*h]; yns=[yns; yn];   % store t's, y's
end
 
plot(tns,yns,'-o')

% Let's redo this for other N's:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:2,

   N=2*N;  h=tmax/N; yn=y0;

   tns=[t0]; yns=[y0]; 

   for i=1:N, 
      yn= yn + h* (a*yn);                 % Euler's method
      tns=[tns; t0+ i*h]; yns=[yns; yn];  % store t's, y's
   end
 
   plot(tns,yns,'-o')
end

title('Euler''s method converges to exact solution as h -> 0 ')
hold off   % release plot