A new algorithm is presented for the numerical analysis of multi-server queueing systems with interarrival and service time distributions of phase type, including both finite and infinite capacity models. The algorithm is based on the iterative solution of balance equations by using successive overrelaxation and aggregation. The main differences with Takahashi's aggregation-disaggregation algorithm [see \textit{Y. Takahashi}, Adv. Appl. Probab. 13, 619-630 (1981; Zbl 0463.60083)] are the use of a dynamically adjusted relaxation factor and the simplified structure due to an adaptation of the disaggregation step. Like Takahashi's method, the algorithm is remarkably robust and the number of iterations required is quite insensitive to the number of states and to the starting point. The performance of the algorithm with a dynamic relaxation factor was found to be superior to that of the conventional aggregation-disaggregation method. Also, in this paper we show how the state probabilities computed by the algorithm can be used to develop a good approximation for the waiting time distribution. Computational experience is given.