Logical Indices outside of Array bounds with multiple iterations of bxsfun and vecnorm.

Hi, I am trying to create a code that will later allow me to automate, using a function .m file, the creation of a cube with 'bubbles within' similar to the shape attached.

When running this code;

clc; clear;
r = [5,3,6]; %Define Vector of radius of bubbles to be created
x = [5,4,-12]; %Define x location of bubbles
y = [8,-12,-1];
z = [-6,6,0];

cubelength = 50; %Length of cube to create

[xg, yg, zg] = meshgrid(-(cubelength/2):1:cubelength/2); %Create mesh of cube
Pcube = [xg(:) yg(:), zg(:)];
bubble = @(x,y,z,r) vecnorm(bsxfun(@minus,Pcavitycube',[x y z]')) > r; %Define what bubble is
Pcavitycube1 = Pcube(bubble(x(:,1),y(:,1),z(:,1),r(:,1)) ,:); %Remove bubble 1 from cube
Pcavitycube2 = Pcavitycube1(bubble(x(:,2),y(:,2),z(:,2),r(:,2)) ,:); %Remove bubble 2 
Pcavitycube3 = Pcavitycube2(bubble(x(:,3),y(:,3),z(:,3),r(:,3)) ,:); %etc

cubeshp = alphaShape(Pcube); %Create Alphashape of mesh of original cube
shp = alphaShape(Pcavitycube3); %Create Alphashape of mesh of cube w/ cavities 
[tri,loc] = alphaTriangulation(shp); %Fill alphashape

model = createpde('thermal','steadystate'); %Create PDE model

[sensor,mshCube] = geometryFromMesh(model,loc',tri'); %Create Geometry

figure, pdegplot(model,'FaceAlpha',0.5); %Plot

I receive the following error message, after completing the 'removal' of mesh in the line Pcavitycube1.

The logical indices in position 1 contain a true value outside of the array bounds.

Error in testingfunc (line 15)
Pcavitycube = Pcavitycube(bubble(x(:,2),y(:,2),z(:,2),r(:,2)) ,:);

I Know this is possible to do using a singular line i.e;

Pcavitycube = Pcube( bubble(5, 8, -6, 5) & bubble(4, -12, 6, 3) & bubble(-12, -1, 0, 6),:);

However I wanted to break it down for use with an unknown number/size/location of bubbles.

Anyhelp decoding the error would be grand.