function y = montemd(a,b,d,n,f)
% Inputs: a & b are d-dim hyperbox bounds; n is sample size; f is a function
% Output: Monte Carlo Approximation for integral of f over hyperbox.
   V = prod(b - a); % Volume
   y = 0; x = zeros (d,1);
   for i = 1:n
      x=a+(b-a).*rand(d,1);
      y = y + f(x);
   end
   y = V*y/n;