function R = romint(f,a,b,n,TOL)
% Implements the Romberg algorithm for approximating an integral.
%
% INPUT: f (function to be integrated)
%        a (left endpoint)
%        b (right endpoint)
%        n (max size of array is n+1)
%        TOL (error tolerance)
%
% OUTPUT: R (Romberg array)
%
R = nan(n+1);
x = zeros(1,2^(n-1));
h = b-a;
x(1:2) = [a,b];
R(1,1) = 0.5*h*sum(f(x(1:2)));
for i = 2:n+1
    h = h/2;
    m = 2^(i-2);
    x(1:m) = linspace(a+h,b-h,m);
    R(i,1) = 0.5*R(i-1,1) + h*sum(f(x(1:m)));
    for j = 2:i
        R(i,j) = R(i,j-1) + (1/(4^(j-1)-1))*(R(i,j-1)-R(i-1,j-1));
    end
    if abs(R(i,i)-R(i-1,i-1)) < TOL
        R = R(1:i,1:i);
        return;
    end
end
fprintf('Method failed to converge in row %d.\n\n',n);