function Volume = MC_sphere(nPoints)

nInside = 0;

for i = 1 : nPoints;
    x = 2 * rand - 1;
    y = 2 * rand - 1;
    z = 2 * rand - 1;
    if ( x^2 + y^2 + z^2 < 1 ) nInside = nInside + 1;
    end;
end;

VolumeCube = 8;
Volume = VolumeCube * nInside / nPoints;
Exact = 4/3*pi;

fprintf('N=%5d  Volume=%.5g  Exact=%.5g  Error=%.5g\n', ...
         nPoints, Volume, Exact, abs(Exact - Volume) );