% This script illustrates the solution of a stiff system with Euler's
% method, with a critical stepsize h = 0.2, in the seamingly OK case
% where the sharp transient does not appear in the solution (due to
% the choice of an initial condition equal to the eigenvector
% corresponding to the small eigenvalue).
%
% If we run Euler long enough we see that roundoff errors end up
% making a nonzero contribution to the term that multiplies the first
% eigenvector... and the numerical solution eventually diverges!  (On
% my Mac this happens for t \approx 11.5)
clear all
close all
f = inline('[-10*y(1)+y(2);-y(2)]','t','y');
t0=0; y0=[1; 9]/sqrt(82); % initial condition = second eigenvector
tmax = 20;
tt=linspace(0,tmax,200);
hvalues = [0.05 0.1 0.2 0.4];
for h=hvalues
    N=tmax/h;
    [t,y] = EulerSystem(t0,y0,f,h,N);
    figure;
    % top plot: first component of solution
    subplot(2,1,1)
    hold on
    xlabel('t'); ylabel('y_1'); 
    plot(t,y(:,1),'o-');
    plot(tt,exp(-tt)/9,'r');
    legend('Approximate solution','Exact soluton',-1);
    title(sprintf('First component of approximate solution with h = %g',h)) 
    hold off
    % bottom plot: second component of solution
    subplot(2,1,2)
    hold on
    xlabel('t'); ylabel('y_2'); 
    plot(t,y(:,2),'o-');
    plot(tt,exp(-tt),'r');
    legend('Approximate solution','Exact soluton',-1);
    title(sprintf('Second component of approximate solution with h = %g',h)) 
    hold off
    pause
end