This article focuses on the implementation of the sampling surfaces method for the three-dimensional vibration analysis of layered piezoelectric plates. The sampling surfaces method is based on choosing inside the layers not equally spaced sampling surfaces parallel to the middle surface in order to introduce the displacements and electric potentials of these surfaces as basic plate variables. Such choice of unknowns with the consequent use of Lagrange polynomials in the assumed distributions of displacements, strains, electric potential, and electric field vector through the thicknesses of layers leads to the robust piezoelectric plate formulation. The sampling surfaces are located inside each layer at Chebyshev polynomial nodes that makes it possible to minimize uniformly the error due to the Lagrange interpolation. Therefore, the sampling surfaces formulation can be applied efficiently to the obtaining of benchmark solutions for the free vibration of layered piezoelectric plates, which asymptotically approach the three-dimensional solutions of piezoelectricity as the number of sampling surfaces tends to infinity.