function R = romberg(f, a, b, max_k)

% function R = romberg(f, a, b, max_k)
% Computes the triangular extrapolation table for Romberg integration
% using the composite trapezoid rule, starting with h=b-a
% f:     function name (either a name in quotes, or an inline function)
% a, b:  lower and upper limits of integration
% max_k: the number of extrapolation steps (= number of columns in R, plus one.)
%        max_k=0 will do no extrapolation.
%
% Example: R = romberg('sin',0,pi,1,6)

 R = zeros(max_k+1);
 for j=1:max_k+1
     R(j,1) = trapezoid(f,a,b,2^(j-1));
 end
 for k=2:max_k+1
    for j=k:max_k+1
        R(j,k) = (4^(k-1)*R(j,k-1)-R(j-1,k-1))/(4^(k-1)-1);
    end
 end