We formulate and implement a spectral scheme designed towards solving the Kohn–Sham equations for clusters in this study. This is motivated by the observation that one of the most successful methods for solving the Kohn–Sham equations for periodic systems – the plane wave method – is a spectral method based on eigenfunction expansion. Our spectral solution method allows for efficient calculation of the electronic structure of clusters with high accuracy and systematic convergence properties without the need for any artificial periodicity. The basis functions in our method form a complete orthonormal set and are expressible in terms of spherical harmonics and spherical Bessel functions. We compute the occupied eigenstates of the discretized Kohn–Sham Hamiltonian using a combination of preconditioned block eigensolvers and Chebyshev polynomial filter accelerated subspace iterations. We highlight several algorithmic and computational aspects of our method, including computation of the electrostatics terms and parallelization. We present results from a variety of benchmark calculations employing local and nonlocal pseudopotentials and compare these to the literature. To illustrate the efficacy of our method, we demonstrate computations involving large systems that contain thousands of electrons. Finally, we briefly discuss the use of our method to study clusters with arbitrary point group symmetries