function [r] = romberg(f, a, b, n )
% romberg numerical integration
% f is the function handle of the function to integrate
% a, b are the lower and upper integration limits
% n is the number of rows in the table

% translation of Cheney and Kincaid algorithm
% shift from 0-based to 1-based algorithm
% r(i,j) --> r(i+1,j+1) everywhere

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

for i = 1 : n-1 % loop over rows i+1 to n-1
   
   % recursive trapezoidal rule for column 1 row i+1
   
   h = h / 2;
   sum = 0.0;
   for k = 1 : 2 : 2^i - 1
      sum = sum + f(a + k*h);
   end
   r(i+1,1) = r(i,1)/2.0 + h*sum;
   
   % Now use richardson extrapolation across row i+1
   % to the diagonal
  
   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