An efficient sparse-matrix-based numerical method is constructed to simulate two-dimensional nuclear magnetic resonance spectra of many-spin systems including the effects of chemical exchange and/or relaxation. The method employs efficient numerical time propagation requiring O(2(2n)) operations in the case of an n spin 1/2 system. Pulses are treated with a fast implementation algorithm achieving O(n2(2n)) scaling (case of n spins 1/2). The method is tested in simulations of double-quantum-filter correlation spectroscopy and exchange spectroscopy experiments on five- and seven-spin systems with two sites. Observed scaling is consistent with the analytic predictions.