A method of reduction of system order is introduced, in which a lower-order transfer function/transfer matrix is obtained by minimizing an error function constructed from the time responses of the system and the reduced model. The denominator of the reduced model can be selected either by retaining dominant poles or by the Routh method. The numerator of the reduced model is obtained by minimizing the step response error along with a steady-state constraint. The error function is converted into the frequency domain and the minima operation is carried out in this domain. Since the numerator polynomial selection is independent of the denominator, the method is very general and can be applied for MIMO systems as well as for SISO systems. A numerical example of a multivariable system illustrates the procedure. © 1988 Taylor & Francis Group, LLC.