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.