A finite-element method is presented for determining the vibrational characteristics of a circular cylinder composed of bonded piezoelectric layers. Finite-element modeling occurs in the radial direction only using quadratic polynomials and the variationally derived partial differential equations are functions of the hoop and axial coordinates (θ, z) and time t. Using solution form Q exp {i(ξz + nθ + ωt)}, with Q as the nodal amplitudes, leads to an algebraic eigensystem where any one of the three parameters (n, ξ, ω), the circumferential or axial wave number or natural frequency, can act as the eigenvalue. Integer values always are assigned to n, leaving two possible eigenvalue problems. With ω as the eigenvalue and real values assigned to ξ, the solutions represent propagating waves or harmonic standing vibrations in an infinite cylinder. When ξ is the eigenvalue and real values assigned to ω, this eigensystem admits both real and complex eigendata. Real ξ’s represent propagating waves or harmonic standing vibrations as noted before. Complex conjugate pairs of ξ’s describe end vibrations, which arise when an incident wave impinges upon a free end of a cylindrical bar. They are standing waves whose amplitudes decay sinusoidally or exponentially from the free end into the interior. Two examples are given to illustrate the method of analysis, viz., a solid piezoelectric cylinder of PZT-4 ceramic material and a two-layer cylinder of PZT-4 covering an isotropic material.