function I = adaptcc(f, a, b, tol)
% Calculate integral of f from a to b with max error of tol (default 1e-6).
if ( nargin < 4 ), tol = 1e-6; end % Default tolerance
[I, C07] = cc(f, a, b);
E = (200*abs(I - C07)).^1.5;       % What to use here?
if ( abs(E) > tol )                % Fail. Subdivide
    c = (a + b) / 2;
    I = adaptcc(f, a, c, tol/2) + ...
        adaptcc(f, c, b, tol/2);
end
end

function [C13, C07] = cc(f, a, b)
[x, c13] = chebpts(13, [a,b]);
[~, c07] = chebpts(7, [a,b]);
ff = f(x);
C13 = c13*ff;
C07 = c07*ff(1:2:end);
end





