Finite difference solutions of the wave equation can describe the propagation of waves in heterogeneous media. High-order finite-difference approaches have been used extensively for long-range propagation for seismological applications. For ultrasound imaging, therapy, and neuromodulation, the description of ultrasound propagation must also include attenuation, dispersion, and nonlinearity. Describing these physical phenomena with high-order, stencils has been a persistent challenge. Here we develop a formulation of the nonlinear, attenuating, dispersing, wave equation in system form that is ideally suited for a wide variety of high-order finite-difference stencils. Relaxation functions are modeled in the wave equation by stretching the velocity and pressure variables into a new relaxed coordinate system that describes attenuation, dispersion, and perfectly matched layers using a single formalism. The final system of equations, composed only of first-order derivatives, allows for the implementation of sophisticated stencils. Then, a high-order finite-difference solution of these equations is presented. Solutions of wave propagation are demonstrated for ultrasound imaging applications for deep human abdominal imaging and compared to experimental measurements showing that the numerical method accurately captures the speckle statistics, aberration strength, reverberation strength, and coherence properties. Applications are demonstrated for transcranial neuromodulation, array design, super-resolution, and machine learning.