In queuing models, phase-type servers are very useful since they can be used to approximate any service distributions of real-life telecommunication networks. However, the corresponding Markov chains are in most cases very large and hard to solve without the help of an efficient simulation model. A simulation framework, based on perfect sampling, had already been developed to address the evaluation of large queuing models. Perfect sampling enables us to compute unbiased samples of the stationary distribution based on a backward coupling of the Markov chain trajectories, starting from all possible states. We use an optimized algorithm which consists in computing a set of extremal envelopes which bound all trajectories to obtain the coupling. When the chain is monotone, the envelopes are directly obtained by the trajectories coming from the extremal states. Otherwise, envelopes are upper and lower bounds over the whole set of trajectories. In this article we prove that phase-type systems are structurally non monotone even with the simplest phase-type distributions like Erlang, hyper-exponential or Coxian-type distributions. Then we provide a generic method for the computation of the upper and lower envelopes for general finite phase-type services in queuing networks. This allows the perfect sampling of such systems and its efficiency is illustrated through several examples.