A computational technique for evaluating the spectral dyadic Green’s function utilized in the analysis of multilayer structures composed of periodic bi-anisotropic layers is presented. The Green’s function is achieved with the help of a fully vectorial rigorous matrix formulation based on a generalized transmission line (TL) formulation of Maxwell’s equations. This matrix formulation simplifies the mathematics involved in the modeling of the coupled diffraction orders within a periodic bi-anisotropic layer. To examine the soundness of this modeling, the extracted Green’s function is used in an integral equation for the analysis of metallic gratings on homogeneous/inhomogeneous anisotropic and chiral materials. The resulting integral equation is then solved by the Method of Moments (MoM) formulated in terms of sub-domain basis and Galerkin’s test functions. As verification, several examples are analyzed and the results obtained by this method are compared with those available in the literature or obtained by EM-solvers. The efficiency assessment of this method is carried out by considering its convergence rate and cost-time in terms of truncation orders. It is revealed not only the metallic grating with inhomogeneous periodic bi-anisotropic layers can be analyzed by this method but also analysis time, memory, and CPU occupancies are decreased in comparison with commercial EM-solvers.