f   = @(t,x) t.^2 + x - x.^4/10;
sol = ode45( f, [0,5], 0 );
% Compute and plot the relative residual on a lot of points
t = RefineMesh( sol.x, 40 );
[ z, dotz ] = deval(sol,t);
deltat = dotz./f(t,z)-1;
figure(1),semilogy(t,abs(deltat),'k.'),set(gca,'fontsize',24)
xlabel('t'),ylabel('relative residual'),axis([0,5,1.0e-8,1.0])
opts = odeset( 'Reltol', 1.0e-6, 'Abstol', 1.0e-11 );
sol = ode45( f, [0,5], 0, opts );
% Compute and plot the relative residual on a lot of points
t = RefineMesh( sol.x, 40 );
[ z, dotz ] = deval(sol,t);
deltat = dotz ./ f(t,z)-1;
figure(2),semilogy( t, abs(deltat), 'k.' ), set(gca,'fontsize',24)
xlabel('t'),ylabel('relative residual'),axis([0,5,1.0e-15,1.0e-3])
opts = odeset( 'Reltol', 1.0e-11, 'Abstol', 1.0e-11 );
sol = ode113( f, [0,5], 1, opts );
% Compute and plot the relative residual on a lot of points
t = RefineMesh( sol.x, 40 );
[ z, dotz ] = deval(sol,t);
deltat = dotz ./ f(t,z)-1;
figure(3),semilogy( t, abs(deltat), 'k.' ), set(gca,'fontsize',16)
xlabel('t'),ylabel('relative residual'),axis([0,5,1.0e-15,1.0e-5])