Real earth media disperse and attenuate propagating mechanical waves. This anelastic behavior can be described well by a viscoelastic model. We have developed a finite-difference simulator to model wave propagation in viscoelastic media. The finite-difference method was chosen in favor of other methods for several reasons. Finite-difference codes are more portable than, for example, pseudospectral codes. Moreover, finite-difference schemes provide a convenient environment in which to define complicated boundaries. A staggered scheme of second-order accuracy in time and fourth-order accuracy in space appears to be optimally efficient. Because of intrinsic dispersion, no fixed grid points per wavelength rule can be given; instead, we present tables, which enable a choice of grid parameters for a given level of accuracy. Since the scheme models energy absorption, natural and efficient absorbing boundaries may be implemented merely by changing the parameters near the grid boundary. The viscoelastic scheme is only marginally more expensive than analogous elastic schemes. The efficient implementation of absorbing boundaries may therefore be a good reason for also using the viscoelastic scheme in purely elastic simulations. We illustrate our method and the importance of accurately modeling anelastic media through 2-D and 3-D examples from shallow marine environments.