function?[sfcount]=acsimpson（fabtol）
%此方法為自適應simpson方法求解定積分
%f為函數，a，b分別為積分區間，tol為精度
%q為積分值，fcount為計算函數值次數

%一下為準備輸入數據
f?=?fcnchk（f‘vectorized‘）;???????%定義內聯函數
h=b-a;?????????????%步長
qu=[a?a+h/2?b];????%三個節點
qu（qu==0）=2^-1024;?%避免計算函數值出現NAN的情況
y=feval（fqu）;?????%計算函數節點處的值
s=y*[1?4?1]‘*h/6;??%Simpson積分
fcount=0;??????????%定義初始計數器
tol=tol*15;????????%調整精度

[sfcount]=subsim（fsquyfcounttol）;?%遞歸函數
fcount=fcount+3;???%計數器加上準備輸入數據的三次計算
end

%以下為遞歸函數
function?[sfcount]=subsim（fsquyfcounttol）
%以下為對qu二分，分別對兩個小區間求出Simpson積分
h=（qu（end）-qu（1））/2;
quadd=[qu（1）+h/2?qu（2）+h/2];
quadd（quadd==0）=2^-1024;?%避免計算函數值出現NAN的情況
yadd=feval（fquadd）;
subs=[y（1）?yadd（1）?