function  app=romb4_3(Fun,a,b,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    function  app=romb4_3(Fun,a,b,n);
%%            ROMBERG INTEGRATION
%%   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
%%     n  - number of intervals to subdivide [a,b]
%%
%%   OUTPUT:
%%    app  - approximate integral
%%     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Step 1
h=b-a;
  R(1,1)=h*(feval(Fun,a)+feval(Fun,b))/2;

%Step 3
for i=2:n;
   p=2^(i-2);
   x=a-.5*h+h*(1:p);
   R(2,1)=(R(1,1)+h*sum(feval(Fun,x)))/2;
   for  j=2:i;
      R(2,j)=(4^(j-1).*R(2,j-1)-R(1,j-1))/(4^(j-1)-1);
   end;
   h=h/2;
   R(1,1:i)=R(2,1:i);
end;
app=R(1,n);
