We develop a numerical method for approximating the surface modes of sphere-like nanoparticles in the quasi-static limit, based on an expansion of (the angular part of) the potentials into spherical harmonics. Comparisons of the results obtained in this manner with exact solutions and with a perturbation ansatz prove that the scheme is accurate if the shape deviations from a sphere are not too large. The method allows fast calculations for large numbers of particles, and thus to obtain statistics for nanoparticles with random shape fluctuations. As an application we present some statistics for the distribution of resonances, polariziabilities, and dipole axes for particles with random perturbations.