Here is the complete working code:
clc
clear
S=input('Number of Simulations: ');
N=input('Number of points: ');
fileID = fopen('Volume.txt','w');
for k = 1:S
[X, Y, Z]=sphere;
a = [1 1 1 1.2; 2.293,1,1,1.8];
%Draw sphere #1
figure('color','white');
s1=surf(X*a(1,4)+a(1,1), Y*a(1,4)+a(1,2), Z*a(1,4)+a(1,3),'FaceColor', [1 1 1],'edgecolor','none','FaceAlpha',0.6);
light
axis equal
set(gca,'Color','k')
hold on
%Draw sphere #2
s2=surf(X*a(2,4)+a(2,1), Y*a(2,4)+a(2,2), Z*a(2,4)+a(2,3),'FaceColor', [0 1 0],'edgecolor','none','FaceAlpha',0.5);
daspect([1 1 1])
view(30,10)
xlabel('x')
ylabel('y')
zlabel('z')
% x1range1 = -0.2;
% x1range2 = 4.093;
% y1range1 = -0.8;
% y1range2 = 2.8;
% z1range1 = -0.8;
% z1range2 = 2.8;
x1range1 = -1.0;
x1range2 = 5.0;
y1range1 = -1.0;
y1range2 = 3.0;
z1range1 = -1.0;
z1range2 = 3.0;
hit1 = 0; %Points inside spheres
hit2 = 0;
for i = 1:N
%Define Box Range and Randomize Points
x = (x1range2-x1range1).*rand(N,1)+x1range1;
y = (y1range2-y1range1).*rand(N,1)+y1range1;
z = (z1range2-z1range1).*rand(N,1)+z1range1;
x1=x-1;
y1=y-1;
z1=z-1;
r1=sqrt(x1.^2+y1.^2+z1.^2);
x2=x-2.293;
y2=y-1;
z2=z-1;
r2=sqrt(x2.^2+y2.^2+z2.^2);
ii = r1<=1.2;
jj = r2<=1.8;
hit1 = sum(ii);
hit2 = sum(jj);
plot3(x(ii), y(ii), z(ii), 'w+');
%plot3(x(~ii), y(~ii), z(~ii), 'r+');
plot3(x(jj), y(jj), z(jj), 'g+');
drawnow
end
SpherePercentVolume = (hit1+hit2)/N;
Volume = ((x1range2-x1range1)*(y1range2-y1range1)*(z1range2-z1range1))*SpherePercentVolume
fprintf(fileID,'%6.3f\r\n',Volume);
end
fclose(fileID);
load Volume.txt
average = mean(Volume);
stdev = std(Volume);
fileID = fopen('AvgVolume.txt','w');
fprintf(fileID,'%6s','Average Volume: ');
fprintf(fileID,'%6.3f\r\n',average);
fprintf(fileID,'%6s','Standard Deviation: ');
fprintf(fileID,'%6.3f',stdev);
