SUMMARYTime-domain simulation is essential for both analysis and design of complex systems. Unfortunately, high model fidelity leads to large system size and bandwidths, often causing excessive computation and memory saturation. In response we develop an efficient scheme for large-order linear time-invariant systems. First, the A matrix is block diagonalized. Then, subsystems of manageable dimensions and bandwidth are formed, allowing multiple sampling rates. Each subsystem is then discretized using a O(n s ) scheme, where n s is the number of states. Subsequently, a sparse matrix O(n s ) discrete-time system solver is employed to compute the history of the state and output. Finally, the response of the original system is obtained by superposition. In practical engineering applications, closing feedback loops and cascading filters can hinder the efficient use of the simulation scheme. Solutions to these problems are addressed in the paper. The simulation scheme, implemented as a MATLAB function fastlsim, is benchmarked against the standard LTI system simulator lsim and is shown to be superior for medium to large systems. The algorithm scales close to O(n 2 s ) for a set of benchmarked systems. Simulation of a high-fidelity model (n s ≈ 2200) of the Space Interferometry Mission spacecraft illustrates real world application of the method.