This paper describes the first reliable Newton algorithm for the sequential computation of the set of dominant poles of scalar and multivariable transfer functions of infinite systems. This dominant pole algorithm incorporates a deflation procedure, which is derived from the partial fraction expansion concept of analytical functions of the complex frequency s and prevents the repeated convergence to previously found poles. The pole residues (scalars or matrices), which are needed in this expansion, are accurately computed by a Legendre-Gauss integral solver scheme for both scalar and multivariable systems. This algorithm is effectively applied to the modal model reduction of multivariable transfer functions for two test systems of considerable complexity and containing many distributed parameter transmission lines.