function s=Self_Adaptive_integral(a,b,tol,M)

%input:  a--下限

%        b--下限

%        tol--the tolerance（容差）

%        m--初始设置的步数

h=(b-a)/M;%步距

s=0;

for i=1:M

    x=a+(i-1)*h;

    y=a+i*h;

    to=abs(simpson_integral(x,y,2)+simpson_integral(x,y,1))/10;

    j=1;

    while(to>=tol)  %循环直到to<tol为止

    j=j+1;

    to=(abs(simpson_integral(x,y,2^j)-simpson_integral(x,y,1)))/10;

    %精度测试式

    end

    s=s+simpson_integral(x,y,2^j);

end

function y=f(x)

%  f--被积函数

y=x^5+sin(x)*x^6*sqrt(x);

function s=simpson_integral(a,b,m)

%input: a--下限

%       b--上限

%       m--步数

h=(b-a)/(2*m);

s1=0;

s2=0;

for i=1:(m-1)

    x=a+2*i*h;

    s1=s1+f(@f,x);

end

for i=1:m

    x=a+(2*i-1)*h;

    s2=s2+f(@f,x);

end

s=h/3*(f(@f,a)+f(@f,b)+2*s1+4*s2);