As a rule, problems of wave propagation in finite media with non-uniform spatial distribution of material properties can only be tackled by numerical models. In addition, the modeling of damage features in a material requires the introduction of locally nonlinear and-more important-non-unique equations of state. Using a multiscale approach, we have implemented a non-linear hysteretic stress-strain relation based on the Preisach-Mayergoyz (PM) model, into a numerical elastodynamic finite integration technique program, which has originally been developed for linearly elastic wave propagation in inhomogeneous media. The simulation results show qualitatively good agreement with data of non-linear resonant bar experiments in homogeneously nonlinear and hysteretic media. When the PM density distribution of hysteretic units at the mesoscopic level is not uniform and/or confined to a finite area in stress-stress space, the response at high amplitude excitation tend to deviate from the quasi-analytical results obtained in the case of a uniform PM-space density. Localized microdamage features in an intact medium can be modeled by conceiving finite zones with pronounced hysteretic stress-strain relations within a "linear" surrounding. Forward calculations reveal a significant influence of the amplitude dependent resonance behavior on the location (edge versus center of a bar), the extend (width of the zone) and the degree (density of hysteretic units) of damage. (C) 2003 Elsevier B.V. All rights reserved.