The presence of compliant features in rock causes nonlinear distortion of sound waves propagating through a sample, even at microstrain levels. When performing resonance experiments, one observes nonlinear behavior in different ways: resonant frequency shifts as function of drive amplitude, nonlinear dissipation including hysteresis, asymmetry of the resonance curves, and rich harmonic spectra. Recently, a new interesting feature has been added to these nonlinearity observations: the existence of a slow dynamics conditioning and recovering characteristic. Experimental evidence will be shown in a companion ASA contribution by TenCate et al. The focus is on the mathematical modeling of all of the nonlinear resonance observations, including the slow dynamics characteristics. The model uses a finite-difference time-domain solution of the one-dimensional wave equation with appropriate boundary conditions, and calculates waveforms along the sample. The key part of the model is the description of a modulus which depends on higher-order elastic constants, hysteresis strength, and its own weighted response over previous times. By integrating the history of the modulus using an exponential weighting function, interesting conditioning and recovering effects can be simulated which agree well with the nonlinear and slow time constant observations in rock. [Work supported by DOE/OBES/UCal.]