How to add new geometry in a loop with PDE toolbox

I simply want to add holes to my geometry in a loop using PDE tool box to find the maximum gradient location and add a hole at this point, then repeat up to 100 loops and see where the holes appear.

Decompose Constructive solid 2-D Geometry (decsg) jumbles the edge labels when more edges are introduced. I have a small section dedicated to finding the boundary values for dirichlet conditions and the rest are all neumann so doesn't matter which order they appear.

I need a way to add a new V_n into my geometry after each loop but , However the decsg requires inputs (gm,sf,ns) which are a mixture of lists and strings, so how can i format this to accomodate a new information in a loop? Thank you in advance for reading my question.

%Add holes V_n into domain R1 at point of max gradient
model = createpde(1);

R1 = [2,6,    1,100,100,50,50,1,   2,2,50,50,100,100]'; %Domain inside boundary
%Voids i manually added but want to autonomize
V1 = [2,4,    50,51,51,50     30,30,31,31]';  % <---Void to be subtracted from domain/Question! 
V1 = [V1; zeros( length(R1)- length(V1),1 )];
V2 = [2,4,    40,41,41,40     40,40,41,41]';
V2 = [V2; zeros(4,1)];
V3 = [2,4,    40,41,41,40     60,60,61,61]';
V3 = [V3; zeros(4,1)];
V4 = [2,4,    49,50,50,49     49,49,50,50]';
V4 = [V4; zeros(4,1)];
voids=[V1, V2, V3,V4];
gm=[R1,voids];
sf='R1-(V1+V2+V3+V4)';
ns=char('R1','V1','V2','V3','V4')';
g = decsg(gm,sf,ns);
geometryFromEdges(model,g);

%Finding the LHS boundary Edge label for general purposes
LHScol=[2;1;1;2;100;0;1]; %LHS colum in g
G=g==LHScol; %find array v in g
sums=sum(G,1); %sum the columns leaving largest value in LHS column
[~,LHS]=max(sums); %find max column 

RHScol=[2;100;100;2;50;1;0]; %RHS colum in g
H=g==RHScol; %find array v in g
sums=sum(H,1); %sum the columns leaving largest value in RHS column
[~,RHS]=max(sums); %find max column 

%Finds all edges that are no Dirichlet
Nedges=(1:length(g));
Nedges(Nedges==LHS)=[];
Nedges(Nedges==RHS)=[];
applyBoundaryCondition(model,'dirichlet','Edge',LHS,'u',1);
applyBoundaryCondition(model,'dirichlet','Edge',RHS,'u',0);
applyBoundaryCondition(model,'neumann','Edge',Nedges,'g',0);

%Coefficients for the general PDE form:
%***   m*U''_t + d*U'_t - grad(c*gradU) + aU = f   ***

specifyCoefficients(model,'m',0,...
                          'd',0,...
                          'c',1,...
                          'a',0,...
                          'f',0);
figure(1)%geometry
title('geometry')
pdegplot(model,'EdgeLabels','on')
xlim([-1.2,1.2])
axis equal

mesh=generateMesh(model,'Hmax',10);

figure(2)%mesh grid
pdemesh(model); 
axis equal

results = solvepde(model);

figure(3)%result plot
pdeplot(model,'XYData',results.NodalSolution)
u=results.NodalSolution;
pdeplot(model,'XYData',u,'ZData',u)

%Find max gradient location
v=(1:1:100);
[X,Y] = meshgrid(v);
querypoints = [X(:),Y(:)]';
[gradx,grady]=evaluateGradient(results,querypoints);
gradx=reshape(gradx,size(X));
grady=reshape(grady,size(Y));
grad=gradx.^2+grady.^2;
[row, col]= find(ismember(grad, max(grad(:))))
peak=[row,col]