function y = monte (a,b,n,f)
%------------------------------------------------------------------------
% Usage:       y = monte (a,b,n,f)
%
% Description: Use the Monte Carlo method to integrate an m-dimensional
%              function f(x).
%
% Inputs:      a = m by 1 vector of lower integration limits
%              b = m by 1 vector of upper integration limits
%                  The domain of integration D is the following
%                  rectangular region:
%                
%                  D = {x: a(k) <= x(k) <= b(k),  1 <= k <= m}
%
%              n = number of random points uses to estimate
%                  integrand (n >= 1)
%              f = string containing name of user-supplied function to 
%                  be integrated. The form of f is:
%
%                  function y = f(x)
%
%                  When f is called with the m by 1 vector x, it must
%                  return the scalar value f(x).
%
% Outputs:     y = estimate of integral. 
%
% Notes:       The error in the estimate of the integral (standard
%              deviation) is roughly proportional to 1/sqrt(n).  
%------------------------------------------------------------------------

% Initialize 

   chkvec (a,1,'monte');
   chkvec (b,2,'monte');  
   n = args (n,1,n,3,'monte');
   m = length (a);
   y = 0;
   x = zeros (m,1);
      
% Compute volume of D and estimate average value of f(x) 
  
   mu = prod(b - a);
   if n > 100
      hwbar = waitbar (0,'Monte Carlo Integration');
   end
   for i = 1 : n
      if n > 100
         waitbar (i/n);
      end
      for j = 1 : m
         x(j) = randu (1,1,a(j),b(j));
      end
      y = y + feval (f,x);
   end
   if n > 100
      close (hwbar);
   end
   y = y*mu/n;
%------------------------------------------------------------------------