龙贝格积分的matlab程序

以下是使用Matlab实现龙贝格积分的程序：

function [I, n] = romberg(f, a, b, tol)
% f: 被积函数
% a, b: 积分区间
% tol: 精度要求
% I: 积分值
% n: 迭代次数

h = b - a;
R(1, 1) = h / 2 * (f(a) + f(b));
n = 1;
for k = 2 : 10
    h = h / 2;
    sum = 0;
    for i = 1 : 2^(k-2)
        x = a + (2 * i - 1) * h;
        sum = sum + f(x);
    end
    R(k, 1) = R(k-1, 1) / 2 + h * sum;
    for j = 2 : k
        R(k, j) = (4^(j-1) * R(k, j-1) - R(k-1, j-1)) / (4^(j-1) - 1);
    end
    if abs(R(k, k) - R(k-1, k-1)) < tol
        break;
    end
    n = n + 1;
end
I = R(k, k);