how to plot packing cirles inside a hexagon?

Hi ! Urgent help needed. Anyone here??? Kindly. I have plotted a bunch of adjacent hexagons within two overlapping hexagons. Now i want to convert these small hexagons into circles in the same pattern. Here is the code i made. I want to convert the smaller hexagons(cyan color) of " Radius=250" (as shown in figure) into circles of same radii. Any simple and easy way or changes in this code?

%Code
Radius=250;
j = 0;
i = 1;
Centre_x(i) = 100;
Centre_y(i) = 100;
R=1000;     % radius of macrocell in meters
v_x = 100+R * cos((0:6)*pi/3); %Vertexes
v_y = 100+R * sin((0:6)*pi/3);
v_x_b = v_x;
v_y_b = v_y;
while Centre_x(i) < R+100  
while Centre_y(i) < ((((3)^0.5)*R)/2)+100

    i = i + 1;
    Centre_x(i) = Centre_x(i-1);
    Centre_y(i) = Centre_y(i-1) + (((3)^0.5)*Radius);
end
j = j + 1;
if mod(j,2) == 0
    i = i + 1;
    Centre_x(i) = Centre_x(i-1) + 3*Radius/2;
    Centre_y(i) = 100;
else
    i = i + 1;
    Centre_x(i) = Centre_x(i-1) + 3*Radius/2;
    Centre_y(i) = ((((3)^0.5)*Radius)/2)+100;
end
end
for i =  1 : size(Centre_x,2)

        Tempx(i,1) = Centre_x(i)- Radius;
        Tempx(i,2) = Centre_x(i) - Radius/2;
        Tempx(i,3) = Centre_x(i) + Radius/2;
        Tempx(i,4) = Centre_x(i) + Radius;
        Tempx(i,5) = Centre_x(i) + Radius/2;
        Tempx(i,6) = Centre_x(i) - Radius/2;
        Tempx(i,7) = Centre_x(i) - Radius;

        Tempy(i,1) = Centre_y(i);
        Tempy(i,2) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,3) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,4) = Centre_y(i);
        Tempy(i,5) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,6) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,7) = Centre_y(i);
end
LL1 = size(Tempx,1); 
j = 0;
for i =  1 : LL1   
     in =  inpolygon(Tempx(i,1),Tempy(i,2),v_x_b,v_y_b);
     if in == 1
         j = j +1        
         V_X_New(j,:) = Tempx(i,:);
         V_Y_New(j,:) = Tempy(i,:);
end
end
j = 1;
i = size(Centre_x,2) + 1;
Centre_x(i) = 100 - 3*Radius/2;
Centre_y(i) = ((((3)^0.5)*Radius)/2)+100;
while Centre_x(i) > -R+100
while Centre_y(i) < (((((3)^0.5)*R)/2))+100 

    i = i + 1;
    Centre_x(i) = Centre_x(i-1);
    Centre_y(i) = Centre_y(i-1) + (((3)^0.5)*Radius);
end
j = j + 1;
if mod(j,2) == 0   
   i = i + 1;
   Centre_x(i) = Centre_x(i-1) - 3*Radius/2;
   Centre_y(i) = 100;
else
    i = i + 1;
    Centre_x(i) = Centre_x(i-1) - 3*Radius/2;
    Centre_y(i) = ((((3)^0.5)*Radius)/2)+100;
end
end
for i =  LL1+1 : size(Centre_x,2)
Tempx(i,1) = Centre_x(i)- Radius;
        Tempx(i,2) = Centre_x(i) - Radius/2;
        Tempx(i,3) = Centre_x(i) + Radius/2;
        Tempx(i,4) = Centre_x(i) + Radius;
        Tempx(i,5) = Centre_x(i) + Radius/2;
        Tempx(i,6) = Centre_x(i) - Radius/2;
        Tempx(i,7) = Centre_x(i) - Radius;

        Tempy(i,1) = Centre_y(i);
        Tempy(i,2) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,3) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,4) = Centre_y(i);
        Tempy(i,5) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,6) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,7) = Centre_y(i);
end
LL2 = size(Tempx,1);
j = size(V_X_New,1); 
for i =  LL1+1 : LL2

     in =  inpolygon(Tempx(i,4),Tempy(i,3),v_x_b,v_y_b);
     if in == 1
         j = j +1;

         V_X_New(j,:) = Tempx(i,:);
         V_Y_New(j,:) = Tempy(i,:);

     end     
end
j = 0;
i = size(Centre_x,2) + 1;
Centre_x(i) = 100;
Centre_y(i) = -((((3)^0.5)*Radius))+100;
while Centre_x(i) > -R+100  
while Centre_y(i) > -((((3)^0.5)*R)/2)+100    
    i = i + 1;
    Centre_x(i) = Centre_x(i-1);
    Centre_y(i) = Centre_y(i-1) - (((3)^0.5)*Radius);
end
j = j + 1;
if mod(j,2) == 0   
   i = i + 1;
   Centre_x(i) = Centre_x(i-1) - 3*Radius/2;
   Centre_y(i) = -((((3)^0.5)*Radius))+100;
else
    i = i + 1;
    Centre_x(i) = Centre_x(i-1) - 3*Radius/2;
    Centre_y(i) = -((((3)^0.5)*Radius)/2)+100;
end
end
for i =  LL2+1 : size(Centre_x,2)          
        Tempx(i,1) = Centre_x(i)- Radius;
        Tempx(i,2) = Centre_x(i) - Radius/2;
        Tempx(i,3) = Centre_x(i) + Radius/2;
        Tempx(i,4) = Centre_x(i) + Radius;
        Tempx(i,5) = Centre_x(i) + Radius/2;
        Tempx(i,6) = Centre_x(i) - Radius/2;
        Tempx(i,7) = Centre_x(i) - Radius;                
        Tempy(i,1) = Centre_y(i);
        Tempy(i,2) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,3) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,4) = Centre_y(i);
        Tempy(i,5) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,6) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,7) = Centre_y(i);
end
LL3 = size(Tempx,1); 
j = size(V_X_New,1); 
for i =  LL2+1 : LL3   
     in =  inpolygon(Tempx(i,4),Tempy(i,5),v_x_b,v_y_b);
     if in == 1
         j = j +1;        
         V_X_New(j,:) = Tempx(i,:);
         V_Y_New(j,:) = Tempy(i,:);
end
end
j = 1;
i = size(Centre_x,2) + 1;
Centre_x(i) = 100 + 3*Radius/2;
Centre_y(i) = -((((3)^0.5)*Radius)/2)+100;
while Centre_x(i) < R+100
while Centre_y(i) > -((((3)^0.5)*R)/2)+100   
    i = i + 1;
    Centre_x(i) = Centre_x(i-1);
    Centre_y(i) = Centre_y(i-1) - (((3)^0.5)*Radius);
end
j = j + 1;
if mod(j,2) == 0
    i = i + 1;
    Centre_x(i) = Centre_x(i-1) + 3*Radius/2;
    Centre_y(i) = -((((3)^0.5)*Radius))+100;
else
    i = i + 1;
    Centre_x(i) = Centre_x(i-1) + 3*Radius/2;
    Centre_y(i) = -((((3)^0.5)*Radius)/2)+100;
end
end
for i =  LL3+1 : length(Centre_x)

        Tempx(i,1) = Centre_x(i)- Radius;
        Tempx(i,2) = Centre_x(i) - Radius/2;
        Tempx(i,3) = Centre_x(i) + Radius/2;
        Tempx(i,4) = Centre_x(i) + Radius;
        Tempx(i,5) = Centre_x(i) + Radius/2;
        Tempx(i,6) = Centre_x(i) - Radius/2;
        Tempx(i,7) = Centre_x(i) - Radius;               
        Tempy(i,1) = Centre_y(i);
        Tempy(i,2) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,3) = Centre_y(i) - ((((3)^0.5)*Radius)/2);
        Tempy(i,4) = Centre_y(i);
        Tempy(i,5) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,6) = Centre_y(i) + ((((3)^0.5)*Radius)/2);
        Tempy(i,7) = Centre_y(i);   
end
LL4 = size(Tempx,1); 
j = size(V_X_New,1); 
for i =  LL3+1 : LL4    
     in =  inpolygon(Tempx(i,1),Tempy(i,6),v_x_b,v_y_b);
     if in == 1
         j = j +1;         
         V_X_New(j,:) = Tempx(i,:);
         V_Y_New(j,:) = Tempy(i,:);
end
end
for i =  1 : size(V_X_New,1)
    plot(V_X_New(i,:),V_Y_New(i,:),'c--');
    x = V_X_New(i,6);
    y = V_Y_New(i,6);
    text1=[num2str(i)]; dx = 0; dy = 90; text(x+dx,y-dy,text1);
end
%end

Here is the figure attached, the result of this code.