Results from elastic-wave simulations in a simple model show that for models characterized by a set of layers with sharp boundaries (discontinuous stiffness tensor), traditional finite-difference methods fail to correctly describe the dynamics of the propagation process. The failure comes from the lack of distinction between model and field variables; the same differential operator is applied to discontinuous (model) and continuous (wavefield) components. This problem is solved with a modified high-order finite-difference modeling scheme (dual-operator method) that uses two distinct operators for evaluating the model and field derivatives, a modified version of Virieux-staggered grid, and Shoenberg-Muir averaging relations. The dual-operator method generates a dynamic response that is more accurate than traditional finite-difference schemes and comparable to Haskell-Thomson propagator-matrix schemes, with the advantage that it can be applied to structurally complex models. The method produces stable, accurate results for both solid and liquid layers.