function s=int(f,a,b)

% This program computes the integral of f(x) from x=a to x=b using an
% adaptive Simpson's rule.  It doesn't use MATLAB's vectorization notation. 

%%% useage (examples at end of file)
%   int(@f,a,b)
%   int(@(x)f,a,b)

% Required arguments:
%  f: the function f(x) is assumed to be defined using a function procedure, 
%     or as an anonymous function (see examples at the end of this file)
%  a: lower endpoint
%  b: upper endpoint

%%% tol and L (below) can be adjusted as needed

%%% tol = error tolerance
tol=1e-6;
%%% value similar to default accuracy used by the MATLAB integral command
% tol=1e-10;

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


if a >= b
	error('You need a < b');
end

nf=1;
m=0.5*(a+b);
fm=f(m);
fa=f(a);  fb=f(b);
F=[a m b fa fm fb (b-a)*(fa+4*fm+fb)/6];
A=[];
Lmax=L;
%%% determine scaling
% xM=linspace(a,b,89);
xM=a+(b-a)*rand(89,1);
for ix=1:89
    ft(ix)=abs(f(xM(ix)));
end
fM=mean(ft);
if fM==0
    fM=1;
end
%%% this is the loop for the levels
for il=2:L
    FF=[];
    %%% find midpoints of left and right subintervals and eval f there
    mLv=0.5*(F(:,1)+F(:,2));
    mRv=0.5*(F(:,2)+F(:,3));
    fmLv=zeros(nf,1);
    fmRv=fmLv;
    for ic=1:nf
        fmLv(ic)=f(mLv(ic));
        fmRv(ic)=f(mRv(ic));
    end
    %%% apply simpson to left and right subintervals and compute errors    
    s2Lv=(F(:,2)-F(:,1)).*(F(:,4)+4*fmLv+F(:,5))/6;
    s2Rv=(F(:,3)-F(:,2)).*(F(:,5)+4*fmRv+F(:,6))/6;
    s2v=s2Lv+s2Rv;    
    Errv=abs(s2v-F(:,7))/15;
    tolAv=fM*tol*(F(:,3)-F(:,1))/(b-a);
    %%% sort errors and decide if new subdivision needed
    %%% A is an array containing completed subintervals and Rhomberg value
    values=[(1:nf)'  Errv(:)];
    values2=sortrows(values,2);
    ind=find(values2(:,2)< tolAv(:),1,'last');
    if isempty(ind)==1
        ind=0;
    end
    if ind>0
        v1=values2(1:ind);
        A=[A; F(v1,1) F(v1,3) (16*s2v(v1)-F(v1,7))/15];
    end
    if ind<nf
        v1=values2(ind+1:nf);
        FF=[F(v1,1) mLv(v1) F(v1,2) F(v1,4) fmLv(v1) F(v1,5) s2Lv(v1)];
        FF=[FF; F(v1,2) mRv(v1) F(v1,3) F(v1,5) fmRv(v1) F(v1,6) s2Rv(v1)]; 
    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
    warning('Maximum number of levels reached')
end
%%% add up values from each subinterval
s=sum(A(:,3));


%%%% Example 1:  f=cos(x)  int(f)=sin(x)
%%%% Using function file
%
% function example1
% a=-7; b=8;
% v=int(@f,a,b)   
% exact=sin(b)-sin(a);
% error=(exact-v)/exact
%
% function y=f(x)
% y=cos(x);
%
%%%% Using anonymous function
% a=-7; b=8;
% v=int(@(x)cos(x),a,b)
% exact=sin(b)-sin(a);
% error=(exact-v)/exact
%
%%%% The error in the above, in absolute value, is less than 1e-10

%%%% Example 2:  f=-60*x^2*sin(x^3)*cos(x^3)^19  int(f)=cos(x^3)^20
%%%% Using function file
%
% function example2
% v=int(@f,0,9.1)                      
% a=0; b=9.1; exact=cos(b^3)^20-cos(a^3)^20;
% error=(exact-v)/exact
%
% function y=f(x)
% y=-60*x^2*sin(x^3)*cos(x^3)^19;
%
%%%% Using anonymous function
% v=int(@(x)-60*x^2*sin(x^3)*cos(x^3)^19,0,9.1)
% a=0; b=9.1; exact=cos(b^3)^20-cos(a^3)^20;
% error=(exact-v)/exact
%
%%%% The error in the above, in absolute value, is less than 2e-6

%%%% Comment about relative error
% The scaling is based on a random sampling, and this has a small affect on
% the value of the relative error.  To illustrate, in Example 1 above, the
% error can vary between 6e-11 and 3e-12.

% What to do if it doesn't work (aside from checking on input errors)
%  1.  Change the value of the error tolerance in line 19. 
%  2.  The number of allowed levels can be increased (see line 24).

%  Background
%  This was written for functions that do not lend themselves to MATLAB's 
%  vector notation.  Consequently, this code takes somewhat more computing 
%  time than those using vectorization (such as MATLAB's integral command).  
%  The algorithm is described in the text "Introduction to  Scientific 
%  Computing and Data Analysis" (Holmes, 2016).  

%  It has been tested using MATLAB, versions R2014b through R2016a(beta)
%  version 1: March 3, 2016
%  version 2: March 9, 2016 (revised to speed code up, and added a small 
%  procedure to handle badly scaled functions)





