% This script illustrates the use of 2nd-order AB method
% for the problem 
% y' = -y^2,   y(0)=1
% on the interval [0,10]
% with exact solution
% y(t) = 1/(1+t)
clear all
clc
close all
scrnsz = get(0, 'Screensize');
f = inline('-y.^2','t','y');
t0=0; y0=1;    % initial condition
tmax = 10;
Nvalues = [50 100 200 400];
for N=Nvalues
    h=tmax/N;
    y1 = 1/(1+h);
    [t,y] = AB2(t0,y0,y1,f,h,N);
    figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
    hold on 
    title(sprintf('Approximate AB2 solution with %d points',N)) 
    xlabel('t')
    ylabel('y')
    plot(t,y,'b','LineWidth',2); 
    plot(t,y,'ro','LineWidth',2); 
    legend('AB2 soluton','AB2 points',-1);
    pause
    title(sprintf('Approximate AB2 and exact solutions with %d points',N)) 
    tt=linspace(0,10,201);
    plot(tt,1./(1+tt),'g','LineWidth',2);  % exact solution
    legend('AB2 soluton','AB2 points', 'Exact Solution',-1);
    hold off
    pause
    figure(length(Nvalues)+1)
    hold on 
    title('Logs of Errors') 
    xlabel('t')
    ylabel('log(|AB2-exact|)')
    err = log(abs(y(2:end)-1./(1+t(2:end))));
    plot(t(2:end),err,'LineWidth',2)
    maxerr = max(abs(y-1./(1+t)));
    if N>Nvalues(1)
        orderestimate = log(oldmax/maxerr)/log(oldh/h)
    end
    oldmax = maxerr;
    oldh = h;
    if N<Nvalues(end)
        pause
    end
end
    