Time-dependent queueing models are important since most of real-life problems are time-dependent. We develop a numerical approximation algorithm for the mean, variance and higher-order moments of the number of entities in the system at time t for the Ph(t)/Ph(t)/s/c queueing model. This model can be thought as a reparameterization to the G(t)/GI(t)/s. Our approach is to partition the state space into known and identifiable structures, such as the M(t)/M(t)/s/c or M(t)/M(t)/1 queueing models. We then use the Polya-Eggenberger distribution to approximate certain unknown probabilities via a two-moment matching algorithm. We describe the necessary steps to validate the approximation and measure the accuracy of the model.