function [sol,v,flag]=simpsonadatt(fun,a,b,toll)
% [sol,v,flag]=simpsonadatt(fun,a,b,toll)
hmin=(b-a)/10000;
sol=0;
c=a;
d=b;
v=[a];
flag=0;
while c~=b 
       if d-c < hmin
           flag=1;
           return
       end
       [solp,err]=estimaerr(fun,c,d);
       if err<toll/2*(d-c)/(b-a)
           sol=sol+solp;
           v=[v d];
           c=d;
           d=b;
       else
           d=(c+d)/2;
       end
end

function [I2,err]=estimaerr(fun,a,b)
h=(b-a)/4;
x=[a:h:b];
fx=feval(fun,x);
I1=(b-a)/6*(fx(1)+4*fx(3)+fx(5));
I2=(b-a)/12*(fx(1)+4*fx(2)+2*fx(3)+4*fx(4)+fx(5));
err=abs(I2-I1)/15;
