A fast full-wave Method of Moments (MoM) solution of the governing integral equation is developed to analyze the problem of electromagnetic(EM) wave radiation and scattering from large finite arrays of printed elements in grounded multilayered media. The full-wave solver developed in this study utilizes a multilayered media dyadic Green's function as the kernel of the governing integral equation. The Green's function is efficiently evaluated via an asymptotic closed-form approximation combined with a fast numerical integration method, which significantly reduces the computational cost required to fill the MoM operator matrix. Various iterative solvers, which include a fast matrix-vector multiplication scheme using the fast Fourier transform (FFT), are used to solve the MoM matrix equation to reduce the storage requirements and solving time, and a discrete Fourier transform (DFT)-based preconditioner is also implemented to accelerate the convergence of iterative solvers. The latter approach is henceforth referred to as the preconditioned iterative (PI)-MoM. Additionally, the method is extended to handle arrays with non-rectangular element truncation boundaries. Furthermore, a hybrid uniform geometrical theory of diffraction (UTD)-MoM method is also developed as an alternative full-wave approach. This hybrid method uses a newly-introduced set of UTD-based global basis functions to drastically reduce the number of unknowns, and can thus eliminate the convergence problem of iterative solvers. The full-wave solver in this work is implemented to predict radiation/scattering from arrays of four types of printed elements, namely printed dipoles, microstrip-line fed patch antennas, probe-fed patch antennas and aperture-coupled patch antennas, in grounded multilayered media. For the modeling of the coaxial of probe-fed patch antennas, a special attachment mode is included to enforce the continuity of the current at the probe-to-patch junction, as well as a magnetic frill generator is used to model the coaxial aperture. A detailed description of the PI-MoM and hybrid UTD-MoM method as well as several numerical results are included in this dissertation to validate the accuracy and utility of these methods.