function T = Romberg(fcn, a, b, itMax)
% Compute the integral of a function
% using trapezoid rule with acceleration
% using Richardson extrapolation
% Best result is given by I(end,end)
% Max iteration level is itMax+1
% K. Ming Leung, 02/10/03

T(1,1) = trapezoid(fcn,a,b,1);
T(2,1) = trapezoid(fcn,a,b,2);
T(2,2) = (4*T(2,1)-T(1,1))/3,
for k = 2:itMax
    n = 2^k;
    T(k+1,1) = trapezoid(fcn,a,b,n);
    for j = 2:k+1
        c = 4^(j-1);
        T(k+1,j) = (c*T(k+1,j-1)-T(k,j-1))/(c-1);
    end
end