龙贝格算法 MATLAB实现

function I = Romberg(f, a, b, epsion)
    I = 0;
    h = b-a;
    k = 0;
    m = 0;
    T = zeros(5);
    %下标转换：T^(k)_0 => T(k+1,1)
    T(1,1) = h/2 * (subs(f,a) + subs(f,b));%即T^(0)_0
    delta = 2*epsion;
    while delta > epsion
        k = k + 1;
        h = h / 2;
        tmpSum = 0; j = a + h;
        while j < b
           tmpSum = tmpSum + subs(f, j);
           j = j + 2*h;
        end
        T(k+1, 1) = 0.5*T(k, 1) + h*tmpSum;
        j = 2;
        while j <= k+1
           T(k+1,j) = (4^(j-1)*T(k+1,j-1) - T(k, j-1)) / (4^(j-1)-1); 
           j = j + 1;
        end
        delta = abs(T(k+1, k+1) - T(k, k));
    end
    I = T(k+1,k+1);
    T
end