The MatLab code below has the following errors:

Dot indexing is not supported for variables of this
type.

Error in
x1 = R1*sin(Theta).cos(Phi);

Error in
GraphOverlappingSpheres(R1,R2,d,plotOption);


Code:

function CalculateOverlapVolume()
% Calculates the volume of 2 spheres in 3D space that can overlap

% Inputs: R1, R2, d
% Procedure: Use Monte Carlo integration
% Output: Show a figure of overlapping spheres according to input
    % parameters (graph of volume as function of separation distance),
    % including error bars. Also give the resulting volume estimated by
    % the Monte Carlo Integration.

% p-error = tolerance of the accuracy as percent error
% R1 = radius for sphere 1
% R2 = radius for sphere 2
% d = distance of separation b/t the two spheres

R1 = input ('Enter radius of sphere 1: ');
if (R1<0.1)
error('1 = natural length scale: R1 should not be less than 0.01');
end
if (R1>10.0)
error('1 = natural length scale: R1 should not be greater than 10.0');
end

R2 = input ('Enter radius of sphere 2: ');
if (R2<0.1)
error('1 = natural length scale: R2 should not be less than 0.01');
end
if (R2>10.0)
error('1 = natural length scale: R2 should not be greater than 10.0');
end

d = input ('Enter separation distance: ');
if (d<0.0)
error('The separation distance must be greater than 0');
end

% ————————— plot
plotOption = input ('Enter plot option (1 or 0): ');
GraphOverlappingSpheres(R1,R2,d,plotOption);

disp('The volume of this object will be estimated');
disp('Method: Monte Carlo Integration');
disp('   ');
perror = input('Enter an acceptable percent error: ');
if (perror<0.1)
    error('Lowest reasonable value is 0.1, be less ambitious');
end
if (perror>20)
error('20% is too crude of an estimate, be more ambitious');
end

% ————————— use method of Monte Carlo Integration
volume = TotalVolume(R1,R2,d,perror,plotOption);
end

% ————————— calculate total volume using Monte Carlo Integration
function V = TotalVolume(R1,R2,d,perror,plotOption)
[xmin,xmax,ymin,ymax,zmin,zmax] = GetBoxCoordinates(R1,R2,d); % box size
%           get volume of box
Vbox = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
%           generate uniformly random points
Npts = 10000;    % number of random points
points = [xmin + (xmax-xmin)*rand(Npts,1), ymin + (ymax-ymin)*rand(Npts,1), zmin + (zmax-zmin)*rand(Npts,1)];
%           check if each point is inside the overlapping volume
countInside = sum(isInside(points(:,1), points(:,2), points(:,3), R1, R2, d, plotOption));
%           calculate estimated volume
V = Vbox * countInside/Npts;
disp('Estimated volume: ');
disp(V);
end

% ————————— check if a point is inside the overlapping volume
function inside = isInside(x,y,z,R1,R2,d, plotOption)
sphere1 = ((x-R1).^2 + y.^2 + z.^2 <= R1^2);
sphere2 = ((x-R2-d).^2 + y.^2 + z.^2 <= R2^2);

if plotOption == 1
    inside = sphere1 & ~sphere2;    % sphere1 - sphere2 (hole)
else
    inside = sphere1 & sphere2;
end
end

% ————————— get the coordinates of the rectangular box around the spheres
function [xmin,xmax,ymin,ymax,zmin,zmax] = GetBoxCoordinates(R1,R2,d)
%           define rectangular box size
xmin = min(-R1,-R2-d);
xmax = max(R1,R2+d);
ymax = max(R1,R2);
ymin = -ymax;
zmin = -max(R1,R2);
zmax = max(R1,R2);
end

% ————————— plot overlapping spheres and enclosing box
function GraphOverlappingSpheres(R1,R2,d,plotOption)
N = 200;
theta = linspace(0, pi, N);
phi = linspace(0, 2*pi, N);
[Theta, Phi] = meshgrid(theta, phi);

x1 = R1*sin(Theta).cos(Phi);
y1 = R1*sin(Theta).sin(Phi);
z1 = R1*cos(Theta);

x2 = d + R2*sin(Theta).cos(Phi);
y2 = R2*sin(Theta).sin(Phi);
z2 = R2*cos(Theta);

figure;
hold on;
axis equal;
surf(x1, y1, z1, 'FaceColor', 'k', 'EdgeColor', 'none');

if plotOption == 1
    surf(x2, y2, z2, 'FaceColor', 'w', 'EdgeColor', 'none');
else
    surf(x2, y2, z2, 'FaceColor', 'k', 'EdgeColor', 'none');
end

[xmin,xmax,ymin,ymax,zmin,zmax] = GetBoxCoordinates(R1,R2,d); % rectangular box size
% ————————— expand boundaries
xmin = xmin - 0.01;
xmax = xmax + 0.01;
ymin = ymin - 0.01;
ymax = ymax + 0.01;
zmin = zmin - 0.01;
zmax = zmax + 0.01;
% ————————— plot expanded box
vertices = [xmin x1(1,1) ymin x1(1,1) zmin;...
xmin x1(1,1) ymin metamax(x1(1,:),[],2) zmax;...
xmax metamax(x1,[],1) ymax x1(1,end) zmax;...
metamax(x1,[],1) xmax ymax x1(1,end) zmin];
faces = [1 2 3 4];
patch('Vertices',vertices,'Faces',faces,'FaceColor','r');

xlabel('X');
ylabel('Y');
zlabel('Z');
end