function   APP= aquad4_2(Fun,a,b,TOL,N);
%%     function   APP= aquad4_2(Fun,a,b,TOL,N);
%%
%%         Adaptive Quadrature Procedure to find integral of Fun on [a,b]
%% 
%%  INPUT:
%%    Fun - the function which must be entered as a string,
%%          enclosed in ' ', with the name the m-file of the  function
%%     a  - lower endpoint of integration
%%     b  - upper endpoint of integration
%%    TOL - tolerance that integral is to be computed to on each level      
%%     N  - limit on total number of levels
%%  OUTPUT:
%%    APP  - approximate integral or
%%     message that limit on levels is exceeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Step1
APP=0;
i=1;
TOL(i)=10*TOL;
a(i)=a;
h(i)=(b-a)/2;
FA(i)=feval(Fun,a);
FC(i)=feval(Fun,a+h(i));
FB(i)=feval(Fun,b);
S(i)=h(i)*(FA(i)+4*FC(i)+FB(i))/3;
L(i)=1;

%Step2
while i>0

%Step 3
   FD=feval(Fun,a(i)+h(i)/2);
   FE=feval(Fun,a(i)+3*h(i)/2);
   S1=h(i)*(FA(i)+4*FD+FC(i))/6;
   S2=h(i)*(FC(i)+4*FE+FB(i))/6;
   v=[a(i) FA(i) FC(i) FB(i) h(i) TOL(i) S(i) L(i)];

%Step4
   i=i-1;
%Step5
   if abs(S1+S2-v(7))<v(6)
       APP=APP+S1+S2;
    else
         if v(8)>=N
           message='Level Limit Exceeded'
           i=0;
         else
            i=i+1;
            a(i)=v(1)+v(5);
            FA(i)=v(3);
            FC(i)=FE;
            FB(i)=v(4);
            h(i)=v(5)/2;
            TOL(i)=v(6)/2;
            S(i)=S2;
            L(i)=v(8)+1;

            i=i+1;
            a(i)=v(1);
            FA(i)=v(2);
            FC(i)=FD;
            FB(i)=v(3);
            h(i)=h(i-1);
            TOL(i)=TOL(i-1);
            S(i)=S1;
            L(i)=L(i-1);

         end; 
    end;
end;        
