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:‘）
????????????