% Initialization step
h = 0.2; h_vect = []; eul_vect = []; rk2_vect = []; rk4_vect = [];

t=0:0.05:3;
x = -1 - t + 2 * exp(t);           	% The exact solution

N = 10;
for i = 1:N
  h_vect = [h_vect h];               	% The new step size vector
  [s y1] = eul ('ypt', 0, 3, 1, h);  	% The approx. solution for Euler method
  [s y2] = rk2 ('ypt', 0, 3, 1, h);  	% The approx. solution for rk2
  [s y4] = rk4 ('ypt', 0, 3, 1, h);  	% The approx. solution for rk4
  plot(t,x, s,y1,'o', s,y2,'+', s,y4,'x') % Plot all
  %plot(t,x, s,y1,'o', s,y2,'+')          % Plot all
  z = -1 - s + 2 * exp(s);           	% The exact solution at the s points
  eulerror = max(abs(z - y1));         	% The biggest error for eul
  rk2error = max(abs(z - y2));         	% The biggest error for rk2
  rk4error = max(abs(z - y4));         	% The biggest error for rk4
  eul_vect = [eul_vect eulerror];    	% The new error vector
  rk2_vect = [rk2_vect rk2error];    	% The new error vector
  rk4_vect = [rk4_vect rk4error];    	% The new error vector
  h_vect;                             	% Display the vector of step sizes
  eul_vect;                           	% Display the vectors of errors
  rk2_vect;
  rk4_vect;
  h = h / 2;                         	% The step size for the next time
end

% See error
loglog(h_vect, eul_vect, '+', h_vect, rk2_vect, '*', ...
    h_vect, rk4_vect, '^', 'LineWidth', 2);  
grid; 
xlabel('Step size'); ylabel('Maximum error'); 
title ('Maximum error vs. step size')
%gtext('eul'); gtext('rk2'); gtext('rk4');
