We present an approach for P-wave modeling in inhomogeneous transversely isotropic media with tilted symmetry axis (TTI media), suitable for anisotropic reverse-time migration. The proposed approach is based on wave equations derived from first principles — the equations of motion and Hooke's law — under the acoustic TI approximation. Consequently, no assumptions are made about the spatial variation of medium parameters. A rotation of the stress and strain tensors to a local coordinate system, aligned with the TI-symmetry axis, makes it possible to benefit from the simple and sparse form of the TI-elastic tensor in that system. The resulting wave equations can be formulated either as a set of five first-order or as a set of two second-order partial differential equations. For the constant-density case, the second-order TTI wave equations involve mixed and nonmixed second-order spatial derivatives with respect to global, nonrotated coordinates. We propose a numerical implementation of these equations using high-order centered finite differences. To minimize modeling artifacts related to the use of centered first-derivative operators, we use discrete second-derivative operators for the nonmixed second-order spatial derivatives and repeated discrete first-derivative operators for the mixed derivatives. Such a combination of finite-difference operators leads to a stable wave propagator, provided that the operators are designed properly. In practice, stability is achieved by slightly weighting down terms that contain mixed derivatives. This has a minor, practically negligible, effect on the kinematics of wave propagation. The stability of the presented scheme in inhomogeneous TTI models with rapidly varying anisotropic symmetry axis direction is demonstrated with numerical examples.