function R = romberg(f, a, b, n)
% Calcula una integral definida en el intervalo [a,b] mediante la
%    tabla de n filas según la fórmula de Romberg.
%
 R = zeros(n,n);
 h = (b-a)./(2.^(0:n-1));
 R(1,1) = (b-a)*(f(a)+f(b))/2;
 for i = 2:n
    subtotal = 0;
    for j=1:2^(i-2)
       subtotal = subtotal + f(a+(2*j-1)*h(i));
    end
    R(i,1) = R(i-1,1)/2+h(i)*subtotal;
    for k=2:i
       R(i,k) = (4^(k-1)*R(i,k-1)-R(i-1,k-1))/(4^(k-1)-1);
    end
 end
end
