  In this paper, we develop a stochastic algorithm based on Euler-Maruyama
scheme to approximate the invariant measure of the limiting multidimensional
diffusion of the $M/Ph/n+M$ queue. Specifically, we prove a non-asymptotic
error bound between the invariant measures of the approximate model from the
algorithm and the limiting diffusion of the queueing model. Our result also
provides an approximation to the steady-state of the diffusion-scaled queueing
processes in the Halfin-Whitt regime given the well established interchange of
limits property. 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 \cite{Gur1} 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 the
$M/Ph/n+M$ queue and its Euler-Maruyama scheme. In particular, the variance of
the CLT of the limiting queue is determined by using Stein's equation and
Malliavin calculus.
