% adaptive Simpson routine
% uses recursive programming

function Sf = adaptiveSimpson(a,b,f,fa,fb,tol,Sf_old)
global pts

c = (a+b)/2;
fc = feval(f,c);
pts = [pts; c];

Sf1 = simp(a,c,f,fa,fc);
Sf2 = simp(c,b,f,fc,fb);
Sf_new = Sf1 + Sf2;

error = abs(Sf_old-Sf_new)/15;
%pause
if error < tol
   Sf = Sf_new;
else
   tol2 = tol/2;
   S1 = adaptiveSimpson(a,c,f,fa,fc,tol2,Sf1);
   S2 = adaptiveSimpson(c,b,f,fc,fb,tol2,Sf2);
   Sf = S1 + S2;
end