Abstract: The multilayer surface integral equations (SIEs) for electromagnetic scattering by infinitely long magnetodielectric cylinders with an arbitrary number of layers are derived and solved by the spectral integral method (SIM). Singularity subtraction for the 2-D Green’s function is used to enhance the computation accuracy and achieve exponential convergence. The final matrix equation after discretization is formed in the Fourier spectral domain rather than the spatial domain, which greatly expedites the SIM solution by accelerating the convolution via the fast fourier transform (FFT) algorithm. A recursive method is proposed to solve the spectral integral equations instead of using an iterative method to lower the computation complexity. Numerical examples for ordinary multilayer cylinders and invisibility cloak cylinders are presented to validate the SIM results by comparing the total fields, scattered fields, and radar cross section (RCS) to analytical solutions or finite-element simulations. They verify that the recursive solution has a complexity of O(MN log N) for an M-layer cylinder with N discrete points on each interface. Meanwhile, the SIM outperforms the analytical method because only the 0th-order and 1st-order special functions (Bessel functions and Hankel functions) are used in the SIM but higher-order functions are necessary for the analytical method to maintain the accuracy.