% This script illustrates the solution of a stiff system with Euler's
% method, with a critical stepsize h = 0.2 (Euler is stable <=> h <
% -2/-|lambda|_max. For this example, |lambda|_max = 10)

clear all
close all
f = inline('[-10*y(1)+y(2);-y(2)]','t','y');
t0=0; y0=[1; 1];    % initial condition
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);
    % y1 = [exp(-h)/9+8*exp(-10*h)/9; exp(-h)]; % y1 from exact solution
    % [t,y] = AB2System(t0,y0,y1,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 + 8*exp(-10*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