function [sol,N,flag]=simpsontoll(fun,a,b,toll)
N=1;
solold=simpson(fun,a,b,N);
flag=1;
for k=1:20
    N=2*N;
    sol=simpson(fun,a,b,N);
    err=abs(sol-solold)/15;
    if err < toll, flag=0; return, end
    solold=sol;
end
