function R = romberg(f, a, b, n)
R(1,1) = .5*(b-a) * (f(a) + f(b));
for j = 2:n
    subtotal = 0;
    hj = (b-a)/2^(j-1);
    for i = 1:2^(j-2)
        subtotal = subtotal + f(a+(2*i-1)*hj);
    end
    R(j,1) = R(j-1,1)/2 + hj*subtotal;
    for k = 2:j
        S = 4^(k-1);
        R(j,k) = (S*R(j, k-1) - R(j-1, k-1)) / (S - 1);
    end
end
