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

global m b k; 
m	=1;
b	=1;
k	= 100;
v0	=1;
x0	=[0.5; 1.5];
t_err	= 1;					%time at which we evaluate the error.
t_end	= 2;					%time until which we integrate the ODE
w0	= sqrt(k/m);

%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
pos2	= x0(2)/w0*sin(w0*t_err) + x0(1)*cos(w0*t_err);	%theroretical position at t=t_err for 2nd system

for i = 1:length(steps)
  tspan		= [0:steps(i):t_end];
  [t,x1]	= my_ode('f1',tspan,v0);
  [t,x2]	= my_ode('f3', tspan, x0);
  idx_err	= find(t==t_err);		%index which corresponds to t=t_err
  err1		= abs(x1(idx_err)-pos1);	
  w0		= sqrt(k/m);
  err2		= abs(x2(idx_err,1)-pos2);	
  newerr	= [err1;err2];			%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');
legend('order1', 'order2');
title('Integration Error for Euler Method');