rombergs method in matlab

main.m
function I = romberg(f, a, b, n)

% Romberg's method for numerical integration
% f: function to integrate
% a, b: limits of integration
% n: number of iterations (at least 2)

h = (b-a)./2.^(0:n-1);

% Compute the first column of the Romberg tableau
R = zeros(n);
R(1,1) = (b-a)*(f(a)+f(b))/2;

% Compute the rest of the tableau
for j = 2:n
    % Compute the terms for this column
    subtotal = 0;
    for i = 1:2^(j-2)
        subtotal = subtotal + f(a + (2*i-1)*h(j));
    end
    R(j,1) = 1/2 * (R(j-1,1) + h(j-1)*subtotal);
    
    % Compute the rest of the column using Richardson extrapolation
    for k = 2:j
        R(j,k) = R(j,k-1) + (R(j,k-1) - R(j-1,k-1))/(4^(k-1) - 1);
    end
end

% Return the most accurate estimate
I = R(n,n);

end
751 chars
33 lines

main.m
f = @(x) sin(x);
a = 0;
b = pi;
n = 4;
I = romberg(f, a, b, n);
64 chars
6 lines