We introduce and successfully test an efficient method to simulate triaxial borehole electromagnetic (EM) induction measurements acquired in axially symmetrical and transversely isotropic (TI) media. The method uses a Fourier series expansion to express the azimuthal dependence of EM fields and the source term whereby the essentially 3D problem collapses to a series of independent 2D problems. Each 2D problem is solved with a semianalytic method that uses normalized Bessel functions and normalized Hankel functions to express the radial dependence of EM fields, thereby improving numerical stability. In addition, use is made of amplitude and slope basis functions to describe the longitudinal dependence of EM fields to avoid grid refinement in the vicinity of horizontal formation boundaries. For validation, we compare the new simulation method to two 1D analytic methods in horizontally and radially layered formations, and to one 3D finite-difference method (3DFD) in multilayered formations that include borehole and invasion zones. Numerical results indicate that the method is accurate in formations with high conductivity contrasts compared to 1D methods and is more than ten times more efficient than the 3DFD method in multilayer formations.