function s = adaptSimp( f, a, b, tol )
  %----------------------------------------------------------------------------
  % J. R. Senning
  % Gordon College
  % February 19, 1999
  % Revised slightly March 2015
  %
  % Usage: s = adaptSimp(@f, a, b, tol)
  %
  % Returns:
  %	s:	The approximate value of the integral of the function f
  %		on the interval [a,b].
  %
  % Parameters:
  %	f:	handle of function to integrate
  %	a:	left endpoint
  %	b:	right endpoint
  %	tol:	desired upper bound on allowable error
  %
  % This function uses an adaptive Simpson's rule to approximate the value of
  % the integral.  Notice that this is not very efficient -- it is recursive
  % and recomputes many function values that have already been computed.  The
  % code in this function is meant to be clear and explain the ideas behind
  % adaptive schemes rather than to be an efficient implementation of one.
  %----------------------------------------------------------------------------

  % Theory says the factor to multiply the tolerance by should be 15, but
  % since that assumes that f'''' is fairly constant, we want to be a bit
  % more conservative...
  tol_factor = 10.0;
  h = 0.5 * ( b - a );
  h2 = 0.5 * h;

  % Compute initial evaluation points corresponding function values
  x0 = a;
  x1 = a + h2;
  x2 = a + h;
  x3 = b - h2;
  x4 = b;
  f0 = f( x0 );
  f1 = f( x1 );
  f2 = f( x2 );
  f3 = f( x3 );
  f4 = f( x4 );

  % Simpson's Rule estimates on single interval and on two half intervals
  s0 = h * ( f0 + 4 * f2 + f4 ) / 3.0;
  s1 = h * ( f0 + 4 * f1 + 2 * f2 + 4 * f3 + f4 ) / 6.0;

  % Either accept weighed average of estimates or split interval in half
  % and sum estimates on each half
  if ( abs( s0 - s1 ) < tol_factor * tol )
    s = s1 + ( s1 - s0 ) / 15.0;
  else
    s = adaptSimp( f, x0, x2, tol/2.0 ) + adaptSimp( f, x2, x4, tol/2.0 );
  end

end
