  In this paper, we develop a stochastic algorithm based on the Euler--Maruyama
scheme to approximate the invariant measure of the limiting multidimensional
diffusion of $G/Ph/n+GI$ queues in the Halfin-Whitt regime. Specifically, we
prove a non-asymptotic error bound between the invariant measures of the
approximate model from the algorithm and the limiting diffusion. To establish
the error bound, we employ the recently developed Stein's method for
multi-dimensional diffusions, in which the regularity of Stein's equation
developed by Gurvich (2014, 2022) plays a crucial role.
  We further prove the central limit theorem (CLT) and the moderate deviation
principle (MDP) for the occupation measures of the limiting diffusion of
$G/Ph/n+GI$ queues and its Euler-Maruyama scheme. In particular, the variances
in the CLT and MDP associated with the limiting diffusion are determined by
Stein's equation and Malliavin calculus, in which properties of a mollified
diffusion and an associated weighted occupation time play a crucial role.
