function [h,S] = compsimptol(a,b,tol)
% Composite Simpson's Rule for function f(x) 
% on [a,b] doubling steps til within tolerance
f = @(x) exp(x);
S0 = 0;
N = 1;
j = 0;
h = (b-a)/2;
S = f(h)*(b-a);
while (abs(S-S0)>tol)
    S0 = S;
    h = (b-a)/(2*N);
    i = 0:N-1;
    xi = a+2*i*h;
    xi1 = a+2*(i+0.5)*h;
    xi2 = a+2*(i+1)*h;
    S = (h/3)*sum(f(xi)+4*f(xi1)+f(xi2));
    j = j + 1;
    N = 2^j;
end

