function monte3
%------------------------------------------------------------------
% Example: Monte Carlo Multidimensional Integration                           
%------------------------------------------------------------------
   clc
   clear
   d = 3; %dimension
   M = 6; % power of 10
   I0 = 1.0; % Exact Answer
   a = -10.*ones(d,1); % Lower Side Vector
   b = 10.*ones(d,1); %Upper Side Vector
   f = inline ('exp(-sum(x.*x)/2) / sqrt(2*pi)^length(x)','x');
   fprintf ('\nMonte Carlo Multidimensional Integration: Two Samples \n');
   fprintf ('\n K    n=10^K        I1        |Error1|        I2        |Error2|\n');
   fprintf ('------------------------------------------------------------------\n');  
   for i = 0 : M
      n = 10^i;
      I1 = montemd(a,b,d,n,f);
      Error1 = 100*(I1-I0)/I0;
      I2 = montemd(a,b,d,n,f);
      Error2 = 100*(I2-I0)/I0;
      fprintf ('%2i %9i %12.6f %12.6f %12.6f %12.6f\n',i,n,I1,Error1,I2,Error2);
   end
   fprintf ('------------------------------------------------------------------\n'); 
%------------------------------------------------------------------
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;
%------------------------------------------------------------------