function [q fnct] = quadadapt(f,a,b,tol,varargin)
% Evaluates definite integral of f(x) from a to b
 global fnct
 if nargin < 4 || isempty(tol), tol = 1.e-6; end
 fnct=0;
 c = (a + b)/2;
 fa = f(a,varargin{:});
 fc = f(c,varargin{:});
 fb = f(b,varargin{:}); fnct=fnct+3;
 q = quadstep(f, a, b, tol, fa, fc, fb, varargin{:});
end

function q = quadstep(f,a,b,tol,fa,fc,fb,varargin)
% Recursive subfunction used by quadadapt.
 global fnct
 h = b - a; c = (a + b)/2;
 fd = f((a+c)/2,varargin{:});
 fe = f((c+b)/2,varargin{:}); fnct=fnct+2;
 q1 = h/6 * (fa + 4*fc + fb);
 q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);
 if abs(q2 - q1) <= tol
    q = q2 + (q2 - q1)/15;
 else
    qa = quadstep(f, a, c, tol, fa, fd, fc, varargin{:});
    qb = quadstep(f, c, b, tol, fc, fe, fb, varargin{:});
    q = qa + qb;
 end
end
