function r = romberg( f, a, b, n )
  % ----------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 17, 1999
  %
  % Usage: r = romberg(@f, a, b, n)
  %
  %    where r        is the Romberg integration array; most accurate
  %                   answer should be at bottom of right-most column,
  %          f        is the handle (precede with @) of the function to
  %                   integrate,
  %          [a, b]   is the interval of integration, and
  %          n        is the number of levels of recursion to use
  %
  % Based on an algorithm in "Numerical Mathematics and Computing"
  % 4th Edition, by Cheney and Kincaid, Brooks-Cole, 1999.
  %
  % The algorithm in the text indexes R(n,m) from 0,0 but Octave starts
  % indices at 1 so that is reflected in the code below.
  % ----------------------------------------------------------------------

  r = zeros( n + 1, n + 1 );
  h = b - a;
  r(1,1) = 0.5 * h * ( f( a ) + f( b ) );

  for i = 1 : n

    % Compute the halved stepsize and use this to sum the function at
    % all the new points (in between the points already computed)
    h = 0.5 * h;
    sum = 0.0;
    for k = 1 : 2 : 2^i - 1
      sum = sum + f( a + k * h );
    end

    % Compute the composite trapezoid rule for the next level of
    % subdivision.  Use Richardson extrapolation to refine these values
    % into a more accurate form.
    r(i+1,1) = 0.5 * r(i,1) + sum * h;
    for j = 1 : i
      r(i+1,j+1) = r(i+1,j) + ( r(i+1,j) - r(i,j) ) / ( 4^j - 1 );
    end

  end

end
