% This script illustrates the solution of a stiff system 
% with critical stepsize h=0.2 by 
% of second-order AB and BDF methods
clear all
clc
close all
scrnsz = get(0, 'Screensize');
f = @(t,y) [-10*y(1)+y(2);-y(2)];
t0=0; y0=[1; 1];    % initial condition
tmax = 10;
tol = 1e-6;
Nvalues = [50 100];
for N=Nvalues
    h=tmax/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);
%    [t,y] = BDF2System(t0,y0,y1,f,h,N,tol);
    figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
    hold on
    xlabel('t')
    ylabel('x')
    plot(t,y(:,1),'b','LineWidth',2); 
    plot(t,y(:,1),'ro','LineWidth',2);
    tt=linspace(0,10,201);
    plot(tt,exp(-tt)/9 + 8*exp(-10*tt)/9,'g','LineWidth',2);  % exact solution
    legend('Approximate solution','Approximate points','Exact soluton',-1);
    title(sprintf('Approximate solution with %d points, x-coordinate',N)) 
    hold off
    pause
    figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
    hold on
    xlabel('t')
    ylabel('y')
    plot(t,y(:,2),'b','LineWidth',2); 
    plot(t,y(:,2),'ro','LineWidth',2); 
    plot(tt,exp(-tt),'g','LineWidth',2);  % exact solution
    legend('Approximate solution','Approximate points','Exact soluton',-1);
    title(sprintf('Approximate solution with %d points, y-coordinate',N)) 
    hold off
    pause
end
    