function [x,y]=trapzoid(x0,y0,x_end,h,jiter,fcn)
%
% function [x,y]=trapzoid(x0,y0,x_end,h,jiter,fcn)
%
% Solve the initial value problem
%    y' = f(x,y),  x0 <= x <= b,  y(x0)=y0
% Use trapezoidal method with a stepsize of h.  The user must 
% supply an m-file to define the derivative f, with some name, 
% say 'deriv.m', and a first line of the form
%   function ans=deriv(x,y)
% jiter is the number of times to iterate the trapezoidal method 
% formula.  A sample call would be
%   [t,z]=trapzoid(t0,z0,b,delta,jiter,'deriv')
%
% Output:
% The routine trapzoid will return two vectors, x and y.
% The vector x will contain the node points
%   x(1)=x0, x(j)=x0+(j-1)*h, j=1,2,...,N
% with 
%   x(N) <= x_end,  x(N)+h > x_end
% The vector y will contain the estimates of the solution Y 
% at the node points in x.
%

% Initialize.
numstep=1;
x(1)=x0; y(1)=y0;
for t=x0:h:x_end-h
  yt=y(numstep);
  numstep=numstep+1;
  x(numstep)=t+h;
%
% Assign initial value for iteration.
%
  if numstep == 2
      yti=yt+h*feval(fcn,t,yt);
  else
      yt0=y(numstep-2);
      yti=yt0+2*h*feval(fcn,t,yt);
  end
% Do iteration.
  for i=1:jiter
      yti=yt+h*(feval(fcn,t,yt)+feval(fcn,t+h,yti))/2;
  end
  y(numstep)=yti;
end
