function eulerdemo(h,tf)

% EULERDEMO demonstrates Euler's method solving the
%           initial value problem
%
%              y' = y  with y(0) = 1
%
%           The exact solution is y(t) = exp(t)
%           
%           The syntax is 
%           eulerdemo(stepaize,finaltime)
%
%           The program can be used without parameters, in which case the
%           default values of the step size is 1, and the default final
%           time is 4.
%
%           The program steps through the numerical approximation to the
%           solution of the equation y' = y with the given step size
%           until it reaches the final time.  At each step it shows the
%           Euler approximation, and illustrates the truncation error
%           made at that step as well as the propagation of that error.
%           The total error made at each step is the sum of this
%           truncation error and the propagation of the 
%           truncation errors made at previous steps.  At the same time
%           the numerical values of these errors are typed to the command
%           window. 
% 
%           To see each step it is necessary to depress the spacebar. 

a=1;
if nargin < 2
  tf=4;
end
if nargin < 1
  h=1;
end
t=linspace(-0.5,tf+0.5,25);
y=exp(a*t);
lw = 2;
ms = 20;
col1 = [1 0 1];
col2 = 'g';
ppph(1)=plot(t,y,'linew',lw);
shg
hold on
ppph(1)=plot(0,1,'.','markersize',ms);
grid
yy=1;
tt=0;
A=1;
while (tt<tf)
  ttnew=min([tt+h,tf]);h=ttnew-tt;
  yynew=yy+a*yy*h;
  ph(1) = plot([tt,ttnew],[yy,yynew],'-','color',col1,'linew',lw);
  ph(2) = plot(ttnew,yynew,'.','markersize',ms,'color',col1);
  ph(3) = plot(ttnew,exp(a*ttnew),'.','markersize',ms);
  shg
  truncerr = A*exp(a*ttnew)-yynew;
  disp(['The truncation error is ',num2str(truncerr)])
  yf=A*exp(tf);
  tt=ttnew;yy=yynew;
  ttt = tt+h:h:tf;
  if tt<tf
    A=yynew*exp(-a*ttnew);
    s=linspace(tt,tf,25);
    z=A*exp(a*s);
    pph(1) = plot(s,z,'--','color',col2,'linew',lw);
%    pph(2) = plot(ttt,A*exp(a*ttt),'.','markersize',ms,'color',col2);
    properr=yf-A*exp(a*tf);
    disp(['This propagates to an error of ',num2str(properr)])
    pause
  end
  
end
disp(['The total error at t = ',num2str(tf), ' is ',num2str(exp(a*tf)-yy)])
hold off

