How can I Associate the users to moving Drones and how to MOVE the Drone toward a specific user in matlab?

How to assign the users (distributed randomly) to the moving drones (also distributed randomly) and can the drone know the positions of associated users to move it toward a specific user ??

% Illustration of drones mobility
function mobilityDrone(NumDrone,ro,center)
datetime('now')
NumDrone= 2; % Number of drones base stations
NumUEs=5;   %Number of users
ro=500; % Raduis of layout [m]
center=[0 0]; % Center of circle for the network layout
v =28.8; % velocity of Drone [km/h]
v = v/3.6;% velocity of Drone [m/sec]
Drone_maxAccel = 2; % Drone maximum acceleration 
tI=1; %[sec] , minimum turning update interval, we will assume
maxtheta = (Drone_maxAccel * tI)/v; % max turning angle for the drone
for G = 3:2:21 % direction candidates , min value of G = 3, G=3,5,7,9,11,13,15,17,19,21
gA = (2*maxtheta)/(G-1); % turning angle step for the drone
end
maxtheta = rad2deg(maxtheta);  %[-14.32; -1.432; 0; 1.432; 14.32];
gA = rad2deg(gA);
TurnAngleOption= [-maxtheta -gA 0 gA maxtheta]; % Drone Turning Angle options
pickangle = TurnAngleOption(randperm(length(TurnAngleOption),1)); %to pick a random number 
%from a vector
% Generate ground Base Station =(gNB)
PosgNB_x = [0 0, 0 0].'; % 
PosgNB_y = [0 0, 0 0].';
plot(PosgNB_x.',PosgNB_y.','k*','MarkerSize', 30);       
hold on
%Generate Drone Base Stations (Drone BS)
theta_Drone=2*pi*(rand(NumDrone,1)); % distributed random number of Base stations
g = 0.5 * ro + 0.5 * ro * rand(NumDrone,1); % let the drones deployed away from the center of 
%circle network layout
PosDrone_x=center(1)+g.*cos(theta_Drone); % Initial positions
PosDrone_y=center(2)+g.*sin(theta_Drone);
PosDrone = [PosDrone_x ,PosDrone_y];
% Generate Users which are associated to the Drones BS and gNB
theta1 = rand(NumUEs, 1) * 2*pi; % distributed random number of Users
r1 = ro * sqrt(rand(NumUEs, 1));
PosUE_x = [r1 .* cos(theta1(:)) + center(:,1)]; 
PosUE_y = [r1 .* sin(theta1(:)) + center(:,2)]; 
PosUE = [r1 .* cos(theta1) + center(1),r1 .* sin(theta1(:)) + center(2)]; 
plot(PosDrone(:,1),PosDrone(:,2),'rp','MarkerSize', 30);
hold on
plot(PosUE(:,1),PosUE(:,2),'r.','MarkerSize', 12);       
hold on
%Associate or assigned these Users to the Drones BS and gNB
distances = hypot(PosDrone_x-PosUE(:,1).', PosDrone_y-PosUE(:,2).');
distances1 = hypot(PosgNB_x-PosUE(:,1).', PosgNB_y-PosUE(:,2).');
[~, assigned_BS] = min(distances, [],1);
[~, assigned_BS1] = min(distances1, [],1);
% Plot these users that associated to the Drones BS and gNB
plot([PosUE(:,1) PosDrone_x(assigned_BS)].', [PosUE(:,2) PosDrone_y(assigned_BS)].', ...
      'color',  [0 0 1]);
plot([PosUE(:,1) PosgNB_x(assigned_BS1)].', [PosUE(:,2) PosgNB_y(assigned_BS1)].', ...
      'color',  [.5 .5 .5]);
grid on
hold on    
axis equal
% Plot the Network layout as Circle
t = linspace(0, 2*pi);
plot(ro * cos(t) + center(1),ro * sin(t) + center(2))
  hfig = figure('Color', 'w');
hdots=plot(PosDrone(:,1),PosDrone(:,2),'bp','MarkerSize', 12);       
hold on
axis equal
    direction = rand(NumDrone, 1) * 2 *pi;
T=10;  % Simulation Time [10 sec]
for t1=0:1:T
    if t1<=1 % to make the drones choose random direction at the initial time
        while ishghandle(hfig)
% Determine new drone's locations
[PosDrone, direction] = step(PosDrone,direction, v);
set(hdots, 'XData',PosDrone(:,1), 'YData' ,PosDrone(:,2))
drawnow % to redraw
       end
    end
    if t1>1
        while ishghandle(hfig)
% Determine new drone's locations
[PosDrone, pickangle] = step(PosDrone, pickangle, v);
set(hdots, 'XData',PosDrone(:,1), 'YData',PosDrone(:,2))
drawnow
        end
    end
end
end
function[PosDroneNew, direction,pickangle] = step(PosDrone, direction,pickangle, v)
v =28.8; % velocity of Drone [km/h]
v = v/3.6;% velocity of Drone [m/sec]
Drone_maxAccel = 2; % in m/sec^2 also we will ass ume [2,4,8], Drone maximum acceleration 
tI=1; %sec , minimum turning update interval, we will assume
maxtheta = (Drone_maxAccel * tI)/v; % max turning angle 
for G = 3:2:21 % direction candidates , min value of G = 3
gA = (2*maxtheta)/(G-1); % turning angle step
end
maxtheta = rad2deg(maxtheta);  %[-14.32; -1.432; 0; 1.432; 14.32];
gA = rad2deg(gA);
TurnAngleOption= [-maxtheta -gA 0 gA maxtheta]; % Drone Turning Angle options
pickangle = TurnAngleOption(randperm(length(TurnAngleOption),1)); %to pick a random number 
%from a vector
% Compute the next position of the drones
    DX = [cos(pickangle(:)) .* v,sin(pickangle(:)) .* v];
    PosDroneNew= PosDrone + DX;
    plot(PosDrone,'r');
   plot(PosDroneNew,'k');
    plot(PosDrone,PosDroneNew);
    % Construct rotation matrix in order to make the drones turn on arc
rotaion = [cosd(pickangle), -sind(pickangle); sind(pickangle), cosd(pickangle)];
%  multiply drone's positions to get the rotated xy
xyRotated = rotaion.*([PosDrone(:,1), PosDroneNew(:,2)]);
xr = xyRotated(:,1);
yr = xyRotated(:,2);
hold on
plot([PosDrone,xr],[PosDroneNew,yr], 'r*-', 'MarkerSize', 30, 'LineWidth', 3)
axis equal; 
end