%function prob_7_5()
% 
% Written by:
% -- 
% John L. Weatherwax                2006-05-29
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; 

% Half the timestep and see if the error decreases by the expected (1/2^4):
hArray = 0.1 * 2.^[0:-1:-4];
err = []; ii=1; 
for h=hArray,
  [tout,yout]=myrk4( @prob_7_5_fn, [0, 2*pi], [ 1 0 ], h );
  % use the initial condition as the exact solution at t=2*pi (the solution is periodic)
  err(ii) = norm(yout(end,:) - yout(1,:)); ii=ii+1; 
end
figure; loglog( hArray, err, 'o' ); grid on; 
xlabel( 'log(stepsize)' ); ylabel( 'log(error)' ); 

% Fit a linear model to this error data:
%
alphas = polyfit( log(hArray(:)), log(err(:)), 1 ) 

disp('Problem 7.5 Part (c)');

h=pi/50;
[tout,yout]=myrk4( @prob_7_5_fn, [0, 2*pi], [ 1 0 ], h );
figure; plot( tout, yout, 'o' ); hold on; grid on; 
fprintf('%s required %d steps, with a global error of %10.6e\n','myrk4',length(tout)-1,norm(yout(end,:) - yout(1,:))); 

opts = odeset('reltol',10^(-6),'refine',1); 
[tout,yout]=ode23( @prob_7_5_fn, [0, 2*pi], [ 1 0 ], opts );
plot( tout, yout, 'x' ); 
fprintf('%s required %d steps, with a global error of %10.6e\n','ode23',length(tout)-1,norm(yout(end,:) - yout(1,:))); 

opts = odeset('reltol',10^(-6),'refine',1); 
[tout,yout]=ode45( @prob_7_5_fn, [0, 2*pi], [ 1 0 ], opts );
plot( tout, yout, 'd' ); 
fprintf('%s required %d steps, with a global error of %10.6e\n','ode45',length(tout)-1,norm(yout(end,:) - yout(1,:))); 

opts = odeset('reltol',10^(-6),'refine',1); 
[tout,yout]=ode113( @prob_7_5_fn, [0, 2*pi], [ 1 0 ], opts );
plot( tout, yout, 'd' ); 
fprintf('%s required %d steps, with a global error of %10.6e\n','ode113',length(tout)-1,norm(yout(end,:) - yout(1,:)));