% Adaptive trapezoid rule
% Usage: r = adaptrap(f,a,b,tol) or [r,x] = adaptrap(f,a,b,tol)
%
% Dependence: trap.m
%
% Input:
% f - Matlab inline function 
% a,b - integration interval
% tol - error tolerance
%
% Output:
% r - computed value of the integral
% x - points used by the algorithm
%
% Examples:
% r=adaptrap(@sin,0,1,0.01);
% r=adaptrap(@myfunc,0,1,1e-4);          here 'myfunc' is any user-defined function in M-file
% r=adaptrap(inline('sin(x)'),0,1,1e-4);
% r=adaptrap(inline('sin(x)-cos(x)'),0,1,1e-4);

function [r,x] = adaptrap(f,a0,b0,tol0)
r = 0;
n = 1;
a(1) = a0;
b(1) = b0;
tol(1) = tol0;
app(1) = trap(f,a,b,1);
x = [a0 b0];
p = 2;
while n > 0
    p = p + 1;
    x(p) = (a(n) + b(n)) / 2;
    c = x(p);
    oldapp = app(n);
    app(n) = trap(f,a(n),c,1);
    app(n+1) = trap(f,c,b(n),1);
    if abs(oldapp - app(n) - app(n+1)) < 3 * tol(n)
        r = r + app(n) + app(n+1);
        n = n - 1;
    else
        b(n+1) = b(n);
        b(n) = c;
        a(n+1) = c;
        tol(n) = tol(n) / 2;
        tol(n+1) = tol(n);
        n = n + 1;
    end
end