function [Q,N] = adaptint(f,a,b,tol,fa,fb)
% [Q,N] = adaptint(f,a,b,tol)
% naive adaptive integration using trapezoid and Simpson rule
% find integral of function f onver interval [a,b] with tolerance tol
% N is the number of function evaluations used

N0 = 0;                         % N0 counts function evaluations
if nargin==4                    % if we are not given fa,fb we have to find them
  fa = f(a); fb = f(b); N0 = N0 + 2;
end
c = (a+b)/2; fc = f(c);  
N0 = N0 + 1;
QT = (b-a)*(fa + fb)/2;         % Trapezoid rule
QS = (b-a)*(fa + 4*fc + fb)/6;  % Simpson rule
if abs(QT-QS) <= tol
  Q = QT;                       % return trapezoid value
  N = N0 ;
else
  [Q1,N1] = adaptint(f,a,c,tol/2,fa,fc);
  [Q2,N2] = adaptint(f,c,b,tol/2,fc,fb);
  Q = Q1 + Q2;
  N = N1 + N2 + N0;
end

