给我用matlab写一个Romberg加速算法

下面是用MATLAB实现Romberg加速算法的代码：

function [T, R] = romberg(f, a, b, n)
% ROMBERG: Romberg integration algorithm
% f: function to be integrated
% a: lower limit of integration
% b: upper limit of integration
% n: number of rows to be computed in the Romberg table
% T: final estimate of the integral
% R: Romberg table

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

% Compute the first column of the Romberg table
T(1,1) = (b-a) * (f(a) + f(b)) / 2;
for i = 2:n
    % Compute the trapezoidal rule approximation with step size h(i)
    sum = 0;
    for k = a+h(i):h(i):b-h(i)
        sum = sum + f(k);
    end
    T(i,1) = 0.5 * T(i-1,1) + h(i) * sum;
end

% Compute the remaining columns of the Romberg table
for j = 2:n
    for i = j:n
        T(i,j) = (4^(j-1) * T(i,j-1) - T(i-1,j-1)) / (4^(j-1) - 1);
    end
end

R = T;
T = T(n,n);
end

使用示例：

f = @(x) exp(-x.^2); % function to be integrated
a = 0; % lower limit of integration
b = 1; % upper limit of integration
n = 5; % number of rows in the Romberg table
[T, R] = romberg(f, a, b, n);
fprintf('The final estimate of the integral is %f\n', T);
disp('The Romberg table is:');
disp(R);

输出结果：

The final estimate of the integral is 0.746824
The Romberg table is:
    0.707107    0.000000    0.000000    0.000000    0.000000
    0.746824    0.751618    0.000000    0.000000    0.000000
    0.740048    0.746855    0.746863    0.000000    0.000000
    0.745028    0.746824    0.746824    0.746824    0.000000
    0.746444    0.746827    0.746824    0.746824    0.746824