Downhole well-logging data acquired by modern electromagnetic (EM) tools enables determining not only conductivity but also conductivity anisotropy of reservoirs. The most popular anisotropic model for EM data interpretation is a horizontally layered model with transverse isotropic conductivity in each layer. Such a model ignores the conductivity change caused by the mud-filtrate invasion that often occurs in a permeable layer. The invasion effect can be so strong that the conductivity derived using the 1D model can be substantially affected. We have developed an efficient forward-modeling approach that includes the invasion in the formation model so that the invasion effect can be properly accounted for in data interpretation and inversion. The approach uses a Fourier series expansion for the electric field in the azimuthal direction to take advantage of the invariance of conductivity in this direction. Each harmonic in the expansion is expressed in terms of numerical eigenmodes in the radial direction and exponential functions in the vertical direction. Physically, the latter describes a set of plane waves propagating upward or downward in the vertical direction. This property allows us to use reflection and transmission matrices to couple EM fields from layer to layer, making it highly efficient to simulate EM logging response because the two matrices are computed only once for all logging points. The approach is best suited for a multilayer and transversely isotropic formation in which each layer can have an arbitrary number of radial discontinuities. Numerical experiments demonstrate that the new approach can accurately model the response of induction and propagation tools in various formations. A speedup of two orders of magnitude is obtained in a multilayer case compared with a previous 2D method using a different hybridization strategy.