A method of reduction is introduced which is an improvement over the conventional Routh approximation procedure. Since some amount of arbitrariness is available in the time moment expressions of the conventional Routh approximant, it can be utilized to minimize the step response error while preserving the stability and moment matching property of the conventional Routh approximant. The method is extended to multivariable systems in the frequency domain.