Residual depth migration is a process of converting an image migrated with one slowness model to an image migrated with another slowness model. The kinematic operators of this process are defined implicitly in a pair of nonlinear equations. This nonlinear system can be efficiently solved by first transforming it into the initial value problem of a pair of partial differential equations, and then solving the partial differential equations with a finite-difference scheme. Once the operators are calculated, the residual depth migration can be accomplished with Kirchhoff-type integrals. The proposed algorithm handles general variable-slowness models, and it is applicable to both pre- and post-stack images. A synthetic example of residual depth migration of a post-stack image shows that the algorithm gives kinematically accurate results.