An efficient and accurate numerical scheme is proposed to solve the incompressible Navier–Stokes equations in a bounded cylinder. The scheme is based on a projection method formulated in primitive variables to maintain the incompressibility constraint, with a second-order semi-implicit scheme for the time integration, and a pseudospectral approximation for the space variables. The Chebyshev-collocation method applied in the radial and axial directions, and the Fourier–Galerkin approximation used in the azimuthal direction lead to a sequence of two-dimensional Helmholtz and Poisson equations for every azimuthal coefficient that are solved by a diagonalization technique. Radial expansions are considered in the diameter of the cell in order to avoid clustering about the axis, and the number of points are selected to ensure that r=0 is not a collocation point. A minimal number of regularity conditions are imposed implicitly at the origin by forcing the proper parity of the Fourier expansions in the radial direction. The method has been tested on analytical solutions and compared with other reliable three-dimensional results. The improvements introduced in the treatment of the spatial discretization reduce significantly the difficulty of implementation of the code, and facilitate the use of high resolutions. Different boundary conditions can also be easily implemented.