The shortest queue problem (SQP) is widely applied for load-balancing procedure in telecommunication systems. Now SQP is one of the most intensively studied and many investigators give full consideration to the SQP research for large-scale queueing systems (LSQS) with a huge number of identical devices because of the integration of 5G/6G networks is carried out at a very high speed. In this work we apply the high-order non-uniform grid numerical scheme for the investigation of the LSQS dynamics with identical single-service devices. We assume that identical services provide by devices with exponentially distributed service time and there is a Poisson incoming flow of requests for LSQS services with finite intensity. LSQS implements the service discipline so that each input request provides a random selection from any m-set devices such device that has the s-th shortest (or equivalently, the (m- s) -th longest) queue size. The evolution dynamics LSQS can be considered using the solutions analysis, where the solutions can be obtained by solving a system of differential equations of infinite degree which can be found with the use of the Markov chains method. We consider the control problem for LSQS, which is formulated in the form of the boundary value problem (BVP) for this system of differential equations with a small parameter. We use the truncation procedure for this singularly perturbed BVP and investigate BVP for the finite order system of differential equations. We apply the high-order non-uniform grid scheme for numerical solving of the truncated BVP. The grid scheme demonstrates good convergence of solutions of the singularly perturbed BVP when a small parameter tend to zero. The results of the LSQS dynamics simulation demonstrate that this LSQS with the control capable of serving a huge number of input requests. © 2023, The Author(s), under exclusive license to Springer Nature Switzerland AG.