Now many researchers are paying great attention to the study of large-scale queueing systems (LSQS) with a lager number of servers because of the development of 5G/6G networks and Internet of Things (IoT) is going at a high speed. All researchers understand that the join-the-shortest-queue rule is widely used balancing mechanisms for such queuing systems. In this paper we use numerical methods for analysis of the evolution of LSQS dynamics with $$n \rightarrow \infty $$ single-services, each with its own exponentially distributed service times of mean $$\bar{t} =1/\mu $$, where $$\mu >0 $$ is a service intensity. We use a Poisson incoming flow of requests with the intensity $$n \lambda >0$$. LSQS implements the service discipline so that each input request provides a randomly selection from any m-set servers such server that has the s-th shortest (or equivalently, the $$(m-s)$$ -th longest) queue size, $$1 \le s \le m$$. The evolution dynamics LSQS can be describe using the function $$u:k^{s,m} (t)$$. The function $$u:k^{s,m} (t)$$ can be found by solving a system of differential equations infinite degree which can be obtained using the Markov chains approach. We formulate the singularly perturbed Cauchy problem for this system of differential equations with a small parameter. We use the truncation procedure for this singularly perturbed Cauchy problem and formulate the finite order system of differential equations. We apply a high-order non-uniform grid scheme for numerical solving of the truncated Cauchy problem. We use different sets of small parameters for time-scaling processes analysis for LSQS. The grid scheme demonstrates good convergence of solutions of the singularly perturbed Cauchy problem when a small parameter tend to zero. The results of the numerical simulation show that this LSQS can hold with a high incoming flow of requests. © 2023, The Author(s), under exclusive license to Springer Nature Switzerland AG.