function [ x, y ] = trapezoid ( f, x_range, y_initial, nstep )
%
%  function [ x, y ] = trapezoid ( f, x_range, y_initial, nstep )
%
%  TRAPEZOID uses NSTEP steps of the trapzezoid PECECE method to 
%  estimate Y, the solution of an ODE, at the equally spaced points X in the 
%  range X_RANGE(1) to X_RANGE(2).  
%
%  One step of the predictor, Euler's forward method,
%  is followed by two steps of the corrector.  
%
%  The name of the derivative function is F. 
%
x(1) = x_range(1);
dx = ( x_range(2) - x_range(1) ) / nstep;
y(:,1) = y_initial;

for i = 1 : nstep

  x(i+1) = x(i) + dx;
%
%  Predict.
%
  y_temp = y(:,i) + dx * feval ( f, x(i), y(:,i) );
%
%  Correct two times.
%
  y_temp = y(:,i) + dx * ( feval ( f, x(i+1), y_temp ) ...
                         + feval ( f, x(i), y(:,i) ) ) / 2.0;

  y_temp = y(:,i) + dx * ( feval ( f, x(i+1), y_temp ) ...
                         + feval ( f, x(i), y(:,i) ) ) / 2.0;

  y(:,i+1) = y_temp;

end