pseudocode for adaptive quadrature

%  Set a small initial h

h = ( b - a ) / 100.0

n = 0
q = 0.0
x0 = a

Loop1 to estimate integral from x0 to x0 + h

 if b <= x0 exit with success

%  Estimate the integral and the error.
%  If the error is small enough, accept the estimate, and advance x0.
%  Otherwise, decrease h and try again.
%
  Loop2 to reduce h if necessary

    if n_max <= n ) exit with error

%  Don't go past b!

    if ( b < x0 + h )
      h = ( b - x0 )
      x1 = x0 + h / 2.0
      x2 = b
    else
      x1 = x0 + h / 2.0
      x2 = x0 + h

%  Compute integral and error estimates using 1 and 2 trapezoids.

    q1 = h * (       f(x0)         +       f(x2) ) / 2.0
    q2 = h * ( 0.5 * f(x0) + f(x1) + 0.5 * f(x2) ) / 2.0
    e1 = abs ( 4.0 * ( q2 - q1 ) / 3.0 )
    e2 = abs (       ( q2 - q1 ) / 3.0 )

%  Decide if h can be increased, or is about right, or needs to be reduced.

    if ( 8.0 * e2 <= tol * h / ( b - a ) )
      h = h * 2.0
      q = q + q2
      x0 = x2
      n = n + 1
      break
    elseif ( e2 <= tol * h / ( b - a ) )
      h = h
      q = q + q2
      x0 = x2
      n = n + 1
      break
    else
      h = h / 2.0
      exit error

    end loop1

  end loop2

end