function simpA

%  Integrate f(x) over [a, b] using adaptive Simpson
%  Include cpu time and number of function evaluations

global ifun

a=0; b=1;

%  n = initial number of subintervals
n=4;

%  error tolerance
tol = 1e-8;

% L = max number of allowed levels
L=25;

%%% ifun = number of function evaluations
ifun=0;

tic
nf=1;
m=0.5*(a+b);
fm=f(m);
fa=f(a);  fb=f(b);
s=simp(a,m,b,fa,fm,fb);
F=[a m b fa fm fb s];
A=[];
Lmax=L;
%  this is the loop for the levels
for il=2:L
    FF=[];
    %%% use simpson on each remaining subinterval
    for is=1:nf
        mL=0.5*(F(is,1)+F(is,2));
        fmL=f(mL);
        mR=0.5*(F(is,2)+F(is,3));
        fmR=f(mR);        
        s2L=simp(F(is,1),mL,F(is,2),F(is,4),fmL,F(is,5));
        s2R=simp(F(is,2),mR,F(is,3),F(is,5),fmR,F(is,6));
        s2=s2L+s2R;
        Err=abs(s2-F(is,7))/15; 
        tolA=tol*(F(is,3)-F(is,1))/(b-a);
        %%% determine if the accuracy is attained
        if Err<tolA
            A=[A; F(is,1) F(is,3) (16*s2-F(is,7))/15];
        else
            FF=[FF; F(is,1) mL F(is,2) F(is,4) fmL F(is,5) s2L];
            FF=[FF; F(is,2) mR F(is,3) F(is,5) fmR F(is,6) s2R];
        end    
    end
    nn=size(FF);
    %%% if no remaining subintervals, then leave the for-loop
    nf=nn(1);
    if nf==0
        Lmax=il;
        break
    end
    F=FF;
end

%%% if max number of levels exceeded, just use what was calcuated
if nf>0
    for is=1:nf
        A=[A; FF(is,1) FF(is,3) FF(is,7)];
    end
end

%%% add up values from each subinterval
ss=size(A);
sum=0;
for is=1:ss(1)
    sum=sum+A(is,3);
end
toc


fprintf('\n')
fprintf(' %12.4e  &  %i \n',sum,ifun)
fprintf('\n')

% number of levels and subintevals used
levels_used=Lmax
subintervals=ss(1)



function g=simp(a,m,b,fa,fm,fb)
g=(b-a)*(fa+4*fm+fb)/6;


function y=f(x)
global ifun
a=sqrt(15);
lambda=pi;
y=((2*a*sinh(x)-2*x)/(1+a^2+x^2+2*a*cosh(x))-1)*sin(lambda*x);
ifun=ifun+1;













