求解多项式方程在曲线[a,b]的积分，同样的方法可以求解一般函数的积分。

function [int] =integRob(xi,a,b)
%a,b为上下界
%xi为多项式系数向量
tk=zeros(1,100);
sk=zeros(1,100);
ck=zeros(1,100);
rk=zeros(1,100);
tk(1)=(b-a)*(0.5*polyval(xi,a)+0.5*polyval(xi,b));
for n=2:20
    sum=0;
    for m=1:2^(n-2)
        sum=sum+polyval(xi,a+(2*m-1)*(b-a)/2^(n-1));
    end
    tk(n)=tk(n-1)/2+sum*(b-a)/2^(n-1);
end
for n=1:20
    sk(n)=tk(n+1)*4/3-tk(n)/3;
end
for n=1:20
    ck(n)=sk(n)*16/15-sk(n)/15;
end
for n=2:20
    rk(n)=ck(n)*64/63-ck(n)/63;
     if (abs(rk(n)-rk(n-1))<=10^(-8))%精度控制
        int=rk(n);
        break;
    end
end
end

验证算法：
选取函数进行积分 ∫ 0 1 1 + x 2 d x \int_{0}^{1}{\sqrt{1+x^2}\mathrm{d} }x ∫01​1+x2 ​dx

f=@(x) sqrt(1+x^2);
int=integRobtest(f,0,1)
int =

   1.147793574520714

与精确解相差小于0.000000001，编写算法求积分可行。