clear
format compact

%%
% We consider the IVP u'=sin[(u+t)^2] over 0 <= t <= 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)';

%%
% Now we perform a convergence study of our two Runge--Kutta
% implementations.
n = 50*2.^(0:5)';
err_eu = 0*n;
err_ieu = 0*n;
for j = 1:length(n)
    [t, u] = eulersys(f, [a,b], u0, n(j));
    err_eu(j) = max(abs(u_exact(t) - u));
    [t, u] = ie2(f, [a,b], u0, n(j));
    err_ieu(j) = max(abs(u_exact(t) - u));
end

%%
%
loglog(n, err_eu, '.-')
hold on 
loglog(n,0.05*(n/n(1)).^(-1),'--')
loglog(n, err_ieu, 'o-')
loglog(n, 0.01*(n/n(1)).^(-2), '--')
xlabel('n')
ylabel('inf-norm error')
title('Convergence of the Runge-Kutta methods')
legend('eu', 'ieu', '1/n', '1/n^2', 'location','southwest')
grid on