In this paper, three-dimensional vibration of rectangular Y-cut crystal plate has been investigated. The three displacement components of plate are expanded in series of Chebyshev polynomial multiplied by the boundary function R which makes expansions satisfy the essential boundary conditions along the edges. Chebyshev polynomial series are chosen as admissible functions for its two distinct advantages: One is that it is a set of complete and orthogonal series in the interval [−1, 1]; the other is that it includes constant and proportional terms. The constant term can easily express the whole rigid displacement of the body. The proportional term can easily reflect the shear force effect along the thickness of a finite plate. The maximum energy function of a plate is expressed in terms of Chebyshev polynomial series. The eigenvalue matrix for natural vibration frequencies is obtained by Ritz method and then solved by computer program. Example of an infinite plate excited by thickness-shear deformation parallel to one edge is solved and verified by exact solutions. Other examples of four clamped edges and four simply supported edges rectangular Y-cut crystal plates are carried out. The trial plates are of three different edge lengths and two different thicknesses. The first twenty frequencies of natural free vibration are compared with those from a finite element method. It also shows that the results from present method and finite element method have a good agreement. Besides, for the advantage of the constant and proportional terms in Chebyshev polynomial series, convergence study demonstrates the rapid rate and high efficiency. The frequencies monotonically decrease and approach certain values with the increase in the number of terms of admissible functions. Three terms are sufficient for the requirement of expansion in the thickness direction. Finally, the free vibration of clamped square Y-cut crystal is investigated. Due to the three-dimensional expansions, any relative displacements in any points of the plate body of modes can be determined very easily by back substitution of the eigenvalues. The deflected shapes of first eight modes show the flexural and thickness extensional modes explicitly.