function[]?=?Romberg（abmMd）
disp（‘請(qǐng)輸入分子多項(xiàng)式a，分母多項(xiàng)式b，積分下限m，積分上限M，計(jì)算精度‘）
f?=?poly2sym（a）/poly2sym（b）??????????????????????%用于給用戶顯示被積函數(shù)的形式?
disp（‘利用梯形公式計(jì)算數(shù)值積分的結(jié)果‘）
kk?=?zeros（）;
kk（11）=1/2*（M-m）/1*（subs（f‘x‘m）+subs（f‘x‘M））;?%先存儲(chǔ)首項(xiàng)?
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;??%通項(xiàng)公式計(jì)算各項(xiàng)值?????????
????kk（i+11）=y?%存儲(chǔ)其他項(xiàng)?
????if（abs（1/3*（kk（i+11）-kk（i1）））<=d）??%判斷梯形公式值是否達(dá)到要求????
????????disp（‘The?result?is:‘）
????????kk（）
????????kk（i+11）%梯形值滿足要求，輸出結(jié)果??????????
????????break;
????else
????????s=（4*kk（i+11）-kk（i1））/（4-1）;%構(gòu)造simpson各項(xiàng)?????????
????????kk（i+12）=s?%存儲(chǔ)????????????
????????if（i+1>=3）
????????????if（i+1>=3?&?abs（1/15*（kk（i+12）-kk（i2）））<=d）
????????????????kk（）
????????????????disp（‘The?result?is:‘）
????????????