A method for approximating high order stable systems via lower order stable models is presented. It is an improvement over previous Routh-Padé methods since the denominator polynomial of the model transfer function, even if based on the elements of a Routh array, contains two free parameters. A systematic procedure is illustrated for evaluating, besides the numerator coefficients, such denominator parameters in order to ensure a good fit, in the Padé sense, and a suitable stability margin.