  Open subwavelength cylindrical resonators of finite height are widely used in
various photonics applications. Circular cylindrical resonators are
particularly important in nanophotonics, since they are relatively easy to
fabricate and can be designed to exhibit different resonance effects. In this
paper, an efficient and robust numerical method is developed for computing
resonant modes of circular cylinders which may have a few layers and may be
embedded in a layered background. The resonant modes are complex-frequency
outgoing solutions of the Maxwell's equations with no sources or incident
waves. The method uses field expansions in one-dimensional (1D) ``vertical''
modes to reduce the original three-dimensional eigenvalue problem to 1D
problems, and uses Chebyshev pseudospectral method to compute the 1D modes and
set up the discretized eigenvalue problem. In addition, a new iterative scheme
is developed so that the 1D nonlinear eigenvalue problems can be reliably
solved. For metallic cylinders, the resonant modes are calculated based on
analytic models for the dielectric functions of metals. The method is validated
by comparisons with existing numerical results, and it is also used to explore
subwavelength dielectric cylinders with high-$Q$ resonances and analyze gold
nanocylinders.
