Many modern electronic/optical devices rely on waves such as solar cells, antennae, radar, and lasers. These devices are mostly built on a patterned layered structure. For optimizing and characterizing these devices, numerical simulations play a crucial role. Two integral equation methods will be presented. Firstly, we developed a robust and fast computational method based on boundary integral equations for the Helmholtz equation in periodically patterned multilayered media. This method uses near- and far-field decomposition to avoid using the quasi-periodic Green’s function. By construction, far-field contribution can be compressed using Schur complement with minimal computational cost. Both the 2-D and 3-D numerical results will be presented. Secondly, two of the main challenges in wave simulation in layered media using layered media Green’s function will be discussed. The layered media Green’s function is represented by Sommerfeld integrals and evaluation of the layered media Green’s function for a given density function is accelerated by adapting the free-space fast multipole method by compressing the Sommerfeld integrals with multipole expansion with transformed basis. In the low-frequency regime, the convergence of multipole and local expansion follows the free-space method results in O(N) algorithm. Additionally, a slow convergence issue of Sommerfeld integral representation of layered media Green’s function when the target and source points are close to each other is overcome with a mathematically equivalent directional representation of Sommerfeld integrals.