自适应求积分

累了，只放代码，各位自行理解

function q = quadadapt(f, a, b, tol, varargin)
% 自适应求积分
if nargin<4||isempty(tol), tol = 1.e-6; end
c = (a+b) / 2;
fa = feval(f, a, varargin{:});
fc = feval(f, c, varargin{:});
fb = feval(f, b, varargin{:});
q = quadstep(f, a, b, tol, fa, fc, fb, varargin{:});
end

function q = quadstep(f, a, b, tol, fa, fc, fb, varargin)
% 递归函数
h = b-a; c = (a+b) / 2;
fd = feval(f, (a+c)/2, varargin{:});
fe = feval(f, (c+b)/2, varargin{:});
q1 = h/6 * (fa + 4*fc + fd);
q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fd);
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