This work introduces a numerical algorithm for solving wave propagation in the presence of an imperfect contact between two isotropic, elastic and heterogeneous media. Non-ideal interfaces of general type can be modeled as boundary discontinuities of the displacement u and its first time derivative (the particle velocity v). The stress field is continuous, and the quantity <&kgr;u+&zgr;v>, where the brackets denote discontinuities across the interface, is equal to the corresponding stress component. The specific stiffness &kgr; introduces frequency dependence and phase changes in the interface response. On the other hand, the specific viscosity &zgr; is related to the energy loss. It is shown here that, in the velocity stress formulation of the wave equation, such a model is described by Maxwell relaxation like functions. I compute the reflection and transmission coefficients in terms of the corresponding incident propagation angle and complex moduli, together with the energy dissipated at the interface. This analysis characterizes the properties of the non-ideal interface. The numerical method is based on a domain decomposition technique that assigns a different mesh to each side of the interface. As stated above, the effects of the interface on wave propagation are modeled through the boundary conditions that require a special boundary treatment based on characteristic variables. The algorithm solves the velocity-stress wave equations and two additional first-order differential equations (in two-dimensional space) in the displacement discontinuity. For each mesh, the spatial derivatives normal to the interface are solved by the Chebyshev method, and the spatial derivatives parallel to the interface are computed with the Fourier method. The algorithm allows general material variability. The modeling is applied to the problems of crack and fracture scattering. ¿ American Geophysical Union 1996