I'm surprised that the sphere is characterized as
sphere = @(X) heaviside(ones(size(X,1),1) - X(:,1).*X(:,2).*X(:,3));
For me it seems that this is the complete cube [-1,1] x [-1,1] x [-1,1] from which the samples are chosen.
In my opinion,
sphere = @(X) (X(:,1).^2 + X(:,2).^2 + X(:,3).^2 <= 1)
so that you get your integral by
function_over_sphere = @(X) (X(:,1).^2 + X(:,2).^2 + X(:,3).^2 <= 1).*(X(:,1).^2.*X(:,2).^2.*X(:,3).^2);
volMonteCarlo(function_over_sphere,-1,1,-1,1,-1,1,100000)
To make one more test, you could try to calculate the volume of the sphere with radius 1 by integrating the constant function 1 over the sphere:
sphere = @(X)(X(:,1).^2 + X(:,2).^2 + X(:,3).^2 <= 1)*1
As you know, the result should approximately be V = 4/3*pi.
The example under
should demonstrate how Monte Carlo Integration works and help to understand the code from above.