Conventional time-domain algorithms that solve the augmented Burgers equation describing nonlinear propagation in a relaxing fluid have a high computational cost associated with discretizing thin shocks. The Burgers equation can be expressed using intrinsic coordinates [Hammerton and Crighton, J. Fluid Mech. 252, 585-599 (1993)], in which waveforms remain single valued beyond the point at which waveform steepening renders them multivalued in physical coordinates. Presented here is a fast simulation approach based on intrinsic coordinates with shocks inserted into the multivalued pressure waveform using the equal area rule from weak shock theory, yielding a single-valued waveform in physical coordinates without requiring fine discretization of the shock. A two-stage approach is developed by first solving the modified Burgers equation in intrinsic coordinates out to distances where shocks are no longer expected to be thin. Beyond that distance, the solution is propagated using a time-domain algorithm in physical coordinates without restrictively small discretization in time. The two-stage approach achieves high accuracy and significantly reduced computational cost when compared to time-domain numerical solutions using exclusively physical coordinates. Computational efficiency and accuracy are demonstrated with simulations of N waves propagating in air and shocks with exponential tails in seawater, each involving two relaxation mechanisms.