The central theme of this paper is to apply generalized block-pulse operational matrices to approximate continuous-time systems. Generalized block-pulse operational matrices and the Routh approximation method are used together to find a low-order transfer function to approximate the original high-order transfer function. The Routh approximation method is used to preserve the stability of the original system by first determining the denominator coefficients of the reduced-order system. Generalized block-pulse operational matrices are then applied to determine numerator coefficients of the reduced-order system by optimally matching the unit step responses of the original and reduced-order systems. This new constrained time-domain matching approach not only yields more satisfactory results than previous methods, but also provides a more straightforward and efficient method for the approximation of continuous-time systems.