function result = monte_carlo_nd ( func, dim_num, a, b, eval_num )

%*****************************************************************************80
%
%% monte_carlo_nd() estimates a multidimensional integral using Monte Carlo.
%
%  Discussion:
%
%    Unlike the other routines, this routine requires the user to specify
%    the number of function evaluations as an INPUT quantity.  
%    No attempt at error estimation is made.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Input:
%
%    function FUNC ( DIM_NUM, X ), the name of the function to be evaluated.
%
%    integer DIM_NUM, the spatial dimension.
%
%    real A(DIM_NUM), B(DIM_NUM), the integration limits.
%
%    integer EVAL_NUM, the number of function evaluations.
%
%  Output:
%
%    real RESULT, the approximate value of the integral.
%
  result = 0.0;
    
  for i = 1 : eval_num

    x = rand ( 1, dim_num );

    result = result + feval ( func, dim_num, x );

  end

  volume = prod ( b(1:dim_num) - a(1:dim_num) );

  result = result * volume / eval_num;

  return
end
