This article deals with 2-D modeling of coupled vibrations of finite piezoelectric resonators. A general solution for all the physical quantities in the Cartesian and cylindrical coordinate systems is deduced from the governing equations by expansion in series summation of trigonometric functions of thickness coordinate and trigonometric or Bessel functions of the lateral one. The essential difference between this model and the earlier ones is that instead of expressing mainly in the thickness coordinate and integration through the thickness, the solutions are expressed in the form of double Fourier series augmented by single Fourier or Fourier-Bessel series, which contributes to better satisfy the mechanical and electrical boundary conditions. The dynamic stiffness matrix of the system is developed. Electrical impedances of a typical piezoelectric parallelepiped under stress-free and symmetrical loading conditions and its frequency spectrum for different width-to-thickness ratios are calculated using our model as well as by the finite element method. A comparison shows an excellent agreement. Finally, theoretical and measured electrical impedances of a piezoelectric parallelepiped and a piezoelectric disk are compared and discussed. The 2-D theoretical model proposed here is shown to be accurate and efficient for coupled vibration analysis of piezoelectric resonators and is applicable for any set of finite dimensions and crystal symmetry.