%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% err_step: Display the integration error as a function of the time step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global b; 
b	=3;
v0	=1;
t_err	= 1;					%time at which we evaluate the error.
t_end	= 2;					%time until which we integrate the ODE

%Here are all the different step sizes for which we'll check the error:
steps	= [5e-4 1e-3 2e-3 5e-3 1e-2 2e-2 5e-2 1e-1 2e-1 5e-1];
err	= [];					%The vector that will contain the errors at each timestep
pos1	= v0*exp(-b*t_err);			%theroretical position at t=t_err for 1st system

for i = 1:length(steps)
  tspan		= [0:steps(i):t_end];
  [t,x]		= my_ode('f1',tspan,v0);
  t1		= find(t==t_err);		%index which corresponds to t=t_err
  newerr = abs(x(t1)-pos1);		%new error value
  err = [err, newerr];				%... that we append to the err vector
end
%A simple  'plot' command will not show a clear linear trend: many
%points will be very close to each other since the timesteps we
%chose are roughtly logarithmically-spaced. Hence, a log-log plot
%if much more befitting here. It works exactly as the plot command:
loglog(steps, err, 'x-'); %'x-' just means we want a to display a line plus a 'x' sign at each data point.			
xlabel('step size');
ylabel('error');
title('Integration Error for Euler Method (order 1 ODE)');