function [H,Eh,S,N] = compsimptolegH(a,b,tol)
% Composite Simpson's Rule for function f(x) 
% on [a,b] doubling steps til within tolerance
f = @(x) exp(x);
ex = exp(4)-1;
S0 = 0;
N = 4;
S = 1;
H = [];
Eh = [];
while (abs(S-S0)>tol)
    S0 = S;
    h = (b-a)/(2*N);
    H = [H,h];
    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));
    eh = abs(S-ex);
    Eh = [Eh,eh];
    N = 2*N;
end
p = polyfit(log(H),log(Eh),1)
nn = length(H);
lnH = log(H);
lnE = log(Eh);
plot([lnH(1),lnH(nn)],[(p(1)*lnH(1)+p(2)),(p(1)*lnH(nn)+p(2))],'b-');
hold on
plot(lnH,lnE,'bo');grid

fontlabs = 'Times New Roman';  
xlabel('$\ln(h)$','FontSize',16,'FontName',fontlabs, ...
    'interpreter','latex');  
ylabel('Absolute Error','FontSize',16,'FontName',fontlabs, ...
    'interpreter','latex'); 
mytitle = '';
title(mytitle,'FontSize',16,'FontName', ...
    'Times New Roman','interpreter','latex'); 
set(gca,'FontSize',16);        

print -depsc egH_gr.eps
end
