function [z, n] = simptol(func,a,b,tol)
% Routine for Simpson's 1/3 rule integration to specified tolerance.
%
% USAGE:  [z, n] = traptol('function',a,b,tol)
%
% input:  function = name of external function file.  This file
%             must be able to receive a vector of x data as input 
%	      and produce a vector of y data as output
%         a, b     = limits of integration
%         tol      = absolute error tolerance in answer
% output: z        = value of integral of function between a and b
%         n        = number of intervals required to obtained desired accuracy


% NOTE: in this routine n refers to the number of segments in the
%   interval from a to b

% first integral using 2 segments
n = 2;
h = (b-a)/n;
x = linspace(a,b,n+1);
y = feval(func,x);
z1 = h*simpson(y);

abserror = Inf;
% start looping with doubled N
while abserror > tol
  n = 2*n;
  h = (b-a)/n;
  x = linspace(a,b,n+1);
  y = feval(func,x);
  z2 = h*simpson(y);
  abserror = abs(z2 - z1);
  z1 = z2; % store this for next loop
end

z = z2;