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);