A product operator matrix is proposed to describe scalar couplings in liquid NMR. Combination of the product operator matrix and non-linear Bloch equations is employed to describe effects of chemical shift, translational diffusion, dipolar field, radiation damping, and relaxation in multiple spin systems with both scalar and dipolar couplings. A new simulation algorithm based on this approach is used to simulate NMR signals from dipolar field effects in the presence of scalar couplings. Several typical coupled spin systems with both intra-molecular scalar couplings and inter-molecular dipolar couplings are simulated. Monte Carlo methods are incorporated into simulations as well to analyze diffusion process in these complicated spin systems. The simulated results of diffusion and relaxation parameters and 2D NMR spectra are coincident with the experimental measurements, and agree with theoretical predictions as well. The simulation algorithm presented herein therefore provides a convenient means for designing pulse sequences and quantifying experimental results in complex coupled spin systems.