% This script illustrates the use of AB methods
% 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
warning off
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;
    [t,yEuler] = Euler(t0,y0,f,h,N);
    y1 = 1/(1+h);      % y1 from exact solution
    %y1 = yEuler(2);     % y1 from Euler
    [t,yAB2] = AB2(t0,y0,y1,f,h,N);
    figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
    hold on 
    title(sprintf('Approximate AB solutions with %d points',N)) 
    xlabel('t')
    ylabel('y')
    plot(t,yEuler,'b-.','LineWidth',2); 
    plot(t,yEuler,'ro','LineWidth',2); 
    plot(t,yAB2,'b','LineWidth',2); 
    plot(t,yAB2,'r+','LineWidth',2); 
    legend('Euler solution','Euler points','AB2 solution','AB2 points',-1);
    pause
    title(sprintf('Approximate AB and exact solutions with %d points',N)) 
    tt=linspace(0,10,201);
    plot(tt,1./(1+tt),'g','LineWidth',2);  % exact solution
    legend('Euler solution','Euler points','AB2 soluton','AB2 points', 'Exact Solution',-1);
    hold off
    pause
    figure(length(Nvalues)+1)
    hold on 
    title('Logs of Errors') 
    xlabel('t')
    ylabel('log(|AB-exact|)')
    errEuler = log(abs(yEuler(2:end)-1./(1+t(2:end))));
    errAB2 = log(abs(yAB2(2:end)-1./(1+t(2:end))));
    plot(t(2:end),errEuler,'b-.','LineWidth',2)
    plot(t(2:end),errAB2,'b','LineWidth',2)
    maxerrEuler = max(abs(yEuler-1./(1+t)));
    maxerrAB2 = max(abs(yAB2-1./(1+t)));
    if N>Nvalues(1)
        orderestimateEuler = log(oldmaxEuler/maxerrEuler)/log(oldh/h)
        orderestimateAB2 = log(oldmaxAB2/maxerrAB2)/log(oldh/h)
    end
    oldmaxEuler = maxerrEuler;
    oldmaxAB2 = maxerrAB2;
    oldh = h;
    if N<Nvalues(end)
        pause
    end
end
warning on
    