function result = quade(fnc,a,b,tol) 
%quade evaluates the numerical integral of the function fnc from a to b
%with relative error tol, using composite 64-point Gaussian integration 
%Default is 10^-6 for tolerance (roughly, 6 correct digits)
%
%Infinite range of integration is allowed, that is, inf and -inf are
%possible values for a and b.
%
%Infinite function values are allowed only if they occur at the endpoints.
%
%Do not forget to use the vector form of operators as .*, ./ and .^
%
%Examples:
% quade(inline('1./x'),1,2)
% quade(inline('1./x'),1,2,10^-14) (This is another way of obtaining log(2)
% up to 14 digits
%f=inline('exp(-x.^2)');quade(f,-inf,inf) (Result should be sqrt(pi)
%Written by: Emre Sermutlu - Halil Tanyer Eyyuboglu
%Ref:

warning off MATLAB:divideByZero

if nargin < 3, error('Missing inputs'); end
if nargin ==3, tol = 10^-6; end 
if isnan(a)==1, error('The limit must be a number or inf');end 
if isnan(b)==1, error('The limit must be a number or inf');end 
if b==a, error('Lower and upper limits cannot be the same'); end

if  isfinite(b)==1 & isfinite(a)==1 
    c=(a+b)/2;
elseif  isfinite(b)==0 & isfinite(a)==1
    c=a+sign(b);
    b=Replace1(fnc,c,b,tol);
elseif  isfinite(b)==1 & isfinite(a)==0
    c=b+sign(a);
    a=Replace1(fnc,a,c,tol);
else    c=0;
    a=Replace1(fnc,a,c,tol);
    b=Replace1(fnc,c,b,tol);
end

if isfinite(feval(fnc,a))==0 
    a=Replace2(fnc,a,c,tol);
end
if isfinite(feval(fnc,b))==0 
    b=Replace2(fnc,c,b,tol);
end

result=Fntint(fnc,a,b,tol);

function F = Fntint(fnc,a,b,tol) 
F=0;
MA=[a b];
%abscissa and weight factors for gaussian integration are given below
k=[0.0243502926634244 0.0729931217877990 0.1214628192961206 0.1696444204239928 ...
        0.2174236437400071 0.2646871622087674 0.3113228719902110 0.3572201583376681 ...
        0.4022701579639916 0.4463660172534641 0.4894031457070530 0.5312794640198945 ...
        0.5718956462026340 0.6111553551723933 0.6489654712546573 0.6852363130542332 ...
        0.7198818501716108 0.7528199072605319 0.7839723589433414 0.8132653151227976 ...
        0.8406292962525804 0.8659993981540928 0.8893154459951141 0.9105221370785028 ...
        0.9295691721319396 0.9464113748584028 0.9610087996520537 0.9733268277899110 ...
        0.9833362538846260 0.9910133714767443 0.9963401167719553 0.9993050417357721];
w=[0.0486909570091397 0.0485754674415034 0.0483447622348030 0.0479993885964583 ...
        0.0475401657148303 0.0469681828162100 0.0462847965813144 0.0454916279274181 ...
        0.0445905581637566 0.0435837245293235 0.0424735151236536 0.0412625632426235 ...
        0.0399537411327203 0.0385501531786156 0.0370551285402400 0.0354722132568824 ...
        0.0338051618371416 0.0320579283548516 0.0302346570724025 0.0283396726142595 ...
        0.0263774697150547 0.0243527025687109 0.0222701738083833 0.0201348231535302 ...
        0.0179517157756973 0.0157260304760247 0.0134630478967186 0.0111681394601311 ...
        0.0088467598263639 0.0065044579689784 0.0041470332605625 0.0017832807216964];
p2=100;%Number of times 64-point rule is used
p1=p2/2;
w8=fliplr(w);w7=[w8 w];w1=repmat(w7,p1,1);w2=repmat(w7,p2,1);
k8=-fliplr(k);k7=[k8 k];k1=repmat(k7,p1,1);k2=repmat(k7,p2,1);
while size(MA,1)~=0
    a=MA(1,1);
    b=MA(1,2);
    MA(1,:)=[];    
    h=(b-a)/(2*p1); %half length of a subinterval
    v8=1:p1;v7=a+(2*v8-1)*h;v1=repmat(v7',1,64); %midpoints of subintervals
    x1=v1+k1*h;
    F1=sum(sum(feval(fnc,x1).*w1))*h; %First Estimate
   
    h=(b-a)/(2*p2); %half length of a subinterval
    v6=1:p2;v5=a+(2*v6-1)*h;v2=repmat(v5',1,64);  %midpoints of subintervals
    x2=v2+k2*h;
    F2=sum(sum(feval(fnc,x2).*w2))*h; %Second Estimate
    
    F8=max([F1 F2]);
    if F8~=0
        err=abs((F2-F1)/F8);
    else
        err=0;
    end

    if err<tol | a+(b-a)/10^5==a 
       F=F+F2;
    else
        d=(b-a)/10;
        a1=a:d:b-d;
        b1=a+d:d:b;
        c1=[a1;b1]';
        MA=[MA;c1];
    end
    
    if size(MA,1)==0
        return
    end
    if size(MA,1)>10000
        error('Divergent integral or too small tolerance')
        return
    end
end

function r1=Replace1(fnc,a,b,tol)
if isfinite(a)==1
    epfin=a;
    epinf=b;
else
    epfin=b;
    epinf=a;
end
tol=tol/10;
mycont=1;
estim=0;
af=epfin;
del=1;
eski=eps;
w1=[41 216 27 272 27 216 41];
w2=[41 216 27 272 27 216 82 216 27 272 27 216 41];
for i=1:300
    del=del*1.25;
    bf=af+sign(epinf)*del;
    h=(bf-af)/6; 
    x=af:h:bf;
    F1=sum(w1.*feval(fnc,x))*h/140;
    h=(bf-af)/12; 
    x=af:h:bf;
    F2=sum(w2.*feval(fnc,x))*h/140;
    F8=max([F1 F2]);
        if F8~=0
            err=abs((F2-F1)/F8);
        else
            err=0;
        end
    if err<0.1
        mytest=F2;
    else
         mytest=Fntint(fnc,af,bf,0.1);
    end
    alter(i)=sign(mytest);
    estim=estim+mytest;
    if  abs(mytest)<tol*abs(estim)&abs(mytest)/abs(eski)<0.95
        mycont=mycont+1;
    end
    eski=mytest;
    if mycont==4
        r1=bf-sign(epinf)*del*2;
        break,end 
    af=bf;
    if i==300
    error('Integral is probably divergent.')
    end
end

function r2=Replace2(fnc,a,b,tol)
if isfinite(feval(fnc,a))==1
    epfin=a;
    asym=b;
else
    epfin=b;
    asym=a;
end
tol=tol/10;
mycont=1;
bf=epfin;
del=(epfin-asym)/2;
estim=0;
eski=eps;
w1=[41 216 27 272 27 216 41];
w2=[41 216 27 272 27 216 82 216 27 272 27 216 41];
for i=1:300 
    af=asym+del;
    h=(bf-af)/6; 
    x=af:h:bf;
    F1=sum(w1.*feval(fnc,x))*abs(h)/140;
    h=(bf-af)/12; 
    x=af:h:bf;
    F2=sum(w2.*feval(fnc,x))*abs(h)/140;
    if F2~=0
         err=abs((F2-F1)/F2);
    elseif F1~=0
         err=abs((F2-F1)/F1);
    else
         err=0;
    end
    if err<0.1
        mytest=F2;
    else
       mytest=Fntint(fnc,af,bf,0.1);
    end
    estim=estim+mytest;
    if  abs(mytest)<tol*abs(estim)&abs(mytest)/abs(eski)<0.92
        mycont=mycont+1;
    end
    eski=mytest;
    if mycont==5 
        r2=asym+del*5;
        break,end  
    del=del/1.5;
    if i==300 | af+(bf-af)/100==af  
    error('Integral is probably divergent.')
    end
    bf=af;
end