clear
format compact

%%
% We consider the IVP u'=\sin[(u+t)^2] over 0\le t \le 4, with u(0)=-1. 
f = @(t,u) sin( (t+u).^2 );
a = 0;  
b = 4;
u0 = -1;

%%
% We use the built-in solver |ode113| to construct an accurate
% approximation to the exact solution.
opt = odeset('abstol',5e-14,'reltol',5e-14); 
uhat = ode113(f,[a,b],u0,opt);
u_exact = @(t) deval(uhat,t)';  % callable function

%%
% Now we perform a convergence study of the AB4 code. 
n = 50*2.^(0:5)';
err2 = 0*n;
err4 = 0*n;
for j = 1:length(n)
    [t,u] = ab2(f,[a,b],u0,n(j));
    err2(j) = max(abs(u_exact(t)-u));
    [t,u] = ab4(f,[a,b],u0,n(j));
    err4(j) = max(abs(u_exact(t)-u));
end

%%
% The method should converge as $O(h^4)$, so a log-log scale is appropriate for the errors. 
loglog(n, err4, 'k.-')
hold on
xlabel('n')
loglog(n, err2, 'k.-')
ylabel('error')
title('Convergence of AB4 and AB2')
loglog(n,0.1*(n/n(1)).^(-4),'--')
loglog(n,0.1*(n/n(1)).^(-2),'--')
legend('AB4','AB2','1/n^4','1/n^2','location','best')
grid on