Here's how I'd do it (if I was in a hurry). First pretend all the spheres are cubes, and find their intersection.
inx = [3363.5, 3363.5, 3346.5, 3366.5, 3357.5];
iny = [1195.5, 1169.5, 1177.5, 1182.5, 1182.5];
inz = [21.5, 32.5, 33.5, 12.5, 12.5];
% Find bounding box for region
xlim = [max(inx - r), min(inx + r)];
ylim = [max(iny - r), min(iny + r)];
zlim = [max(inz - r), min(inz + r)];
Then, generate a large number ( ngrid^3 ) of points inside that box
ngrid = 100;
[X, Y, Z] = ndgrid(linspace(xlim(1), xlim(2), ngrid), ...
linspace(ylim(1), ylim(2), ngrid), ...
linspace(zlim(1), zlim(2), ngrid));
X = X(:); Y = Y(:); Z = Z(:);
Finally, rule them out by checking every point against every sphere. If you're lucky there should be some left.
% Now rule points out
idx = true(size(X)); % All points assumed in to start with
for i = 1:numel(inx)
idx = idx & (X - inx(i)).^2 + (Y - iny(i)).^2 + (Z - inz(i)).^2 <= r^2;
end
X = X(idx); Y = Y(idx); Z = Z(idx);
In your case, the intersection of all spheres is empty.
