龙贝格算法MATLAB程序

function[] = Romberg（abmMd）
disp（‘请输入分子多项式a，分母多项式b，积分下限m，积分上限M，计算精度‘）
f = poly2sym（a）/poly2sym（b）                      %用于给用户显示被积函数的形式 
disp（‘利用梯形公式计算数值积分的结果‘）
kk = zeros（）;
kk（11）=1/2*（M-m）/1*（subs（f‘x‘m）+subs（f‘x‘M））; %先存储首项?
for  i=1:1:2^40
    t=0;
    for  j=0:1:2^（i-1）-1
        v=m+（2*j+1）*（M-m）/（2^i）;
        vv=polyval（av）/polyval（bv）;
        t=t+（M-m）/（2^i）*vv;
    end
    y=1/2*kk（i1）+t;  %通项公式计算各项值?????????
    kk（i+11）=y %存储其他项?
    if（abs（1/3*（kk（i+11）-kk（i1）））<=d）  %判断梯形公式值是否达到要求????
        disp（‘The?result?is:‘）
        kk（）
        kk（i+11）%梯形值满足要求，输出结果??????????
        break;
    else
        s=（4*kk（i+11）-kk（i1））/（4-1）;%构造simpson各项?????????
        kk（i+12）=s %存储????????????
        if（i+1>=3）
            if（i+1>=3 & abs（1/15*（kk（i+12）-kk（i2）））<=d）
                kk（）
                disp（‘The result is:‘）