% adaptive trapezoidal rule
% u is previous estimate of integral of f on interval [a,b]
% tol is desired tolerance
function val=AdaptiveTrapRule(u,a,b,tol)

h=(b-a)/2;
uleft=(f(a)+f((a+b)/2))*h/2;
uright=(f((a+b)/2)+f(b))*h/2;
usplit=uleft+uright;
if abs(u-usplit)<tol || h<1e-12
    val=usplit;
else
    val=AdaptiveTrapRule(uleft,a,a+h,tol/2)+AdaptiveTrapRule(uright,a+h,b,tol/2);
end

    function y=f(x)
        y=sin(x);
    end
end