A sparse-matrix-based numerical method is constructed to simulate NMR spectra of many-spin systems, including the effects of chemical exchange and/or relaxation. The associated computational demands are predicted to scale like O(2(2n)), as the number of spins n increases; This is vastly superior to the inevitable O(2(6n)) scaling of conventional Householder-based methodology. The improved scaling is verified via numerical computations of simple isomerization systems with four to nine spins. The new method is based on splitting the propagator and use of Chebyshev polynomial expansion of the exponential function.