function R = romberg(f, a, b, n)
% Compute integral of f(x) from a to b using Romberg integration.

h = b - a;
R = zeros(n, n);
R(1, 1) = h / 2 * (f(a) + f(b));
for i = 2:n
    x = a + ((1:2^(i - 2)) - 0.5) * h;
    R(i, 1) = (R(i - 1, 1) + h * sum(f(x))) / 2;
    for j = 2:i
        R(i, j) = R(i, j - 1) + ...
            (R(i, j - 1) - R(i - 1, j - 1)) / (4^(j - 1) - 1);
    end
    h = h / 2;
end
