  A time-domain numerical code based on the constitutive relations of nonlinear
acoustics for simulating ultrasound propagation is presented. To model
frequency power law attenuation, such as observed in biological media, multiple
relaxation processes are included and relaxation parameters are fitted to both
exact frequency power law attenuation and empirically measured attenuation of a
variety of tissues that does not fit an exact power law. A computational
technique based on artificial relaxation is included to correct the
non-negligible numerical dispersion of the numerical method and to improve
stability when shock waves are present. This technique avoids the use of high
order finite difference schemes, leading to fast calculations. The numerical
code is especially suitable to study high intensity and focused axisymmetric
acoustic beams in tissue-like medium, as it is based on the full constitutive
relations that overcomes the limitations of the parabolic approximations, while
some specific effects not contemplated by the Westervelt equation can be also
studied. The accuracy of the method is discussed by comparing the proposed
simulation solutions to one-dimensional analytical ones, to $k$-space numerical
solutions and also to experimental data from a focused beam propagating in a
frequency power law attenuation media.
