% the program takes the intial data y_0 at the initial time x_0 and 
% computes up to time x_final using N timesteps of size dx.  It creates
% two vectors, x and y, which you can then plot against each other.
% The ODE you're trying to solve is dy/dx = f(x,y) where the vector
% field is contained in the subroutine f.m
%
% To call the routine, use the command:
%
% [x,y] = richardson_extrap_euler_step_3(y_0,x_0,x_final,N,@f)

function [x,y] = richardson_extrap_euler_step_3(y_0,x_0,x_final,N,f)

% step-size:
dx = (x_final-x_0)/N;

% put initial data into vector
x(1) = x_0;
y(:,1) = y_0;
for ii=1:N
  F1 = feval(f,x(ii),y(:,ii));
  % take one big step to reach y_big
  y_big = y(:,ii) + dx*F1;
  % take two small steps to reach y_med
  y_temp = y(:,ii) + dx/2*F1;
  F = feval(f,x(ii)+dx/2,y_temp);
  y_med = y_temp + dx/2*F;
  % take four smaller yet steps to reach y_small
  % y at x + dx/4
  y_temp = y(:,ii) + dx/4*F1;
  F = feval(f,x(ii)+dx/4,y_temp);
  % y at x + 2*(dx/4)
  y_temp = y_temp + dx/4*F;
  F = feval(f,x(ii)+2*(dx/4),y_temp);
  % y at x + 3*(dx/4)
  y_temp = y_temp + dx/4*F;
  F = feval(f,x(ii)+3*(dx/4),y_temp);
  y_small = y_temp + dx/4*F;

  % use Richardson extrapolation
  rich_big = 2*y_med - y_big;
  rich_small = 2*y_small - y_med;

  y(:,ii+1) = (4/3)*rich_small - (1/3)*rich_big;
  x(ii+1) = x_0 + ii*dx;
end