We present a numerically improved multipole formulation for the calculation of resonances of multiple disks located at arbitrary positions in a 2-d plane, and suitable for the accurate computation of the resonances of large numbers of disks and of high-wavenumber eigenstates. Using a simple reformulation of the field expansions and boundary conditions, we are able to transform the multipole formalism into a linear eigenvalue problem, for which fast and accurate methods are available. Observing that the motion of the eigenvalues in the complex plane is analytic with respect to a two parameter family, we present a numerical algorithm to compute a range of multiple-disk resonances and field distributions using only two diagonalizations. This method can be applied to photonic molecules, photonic crystals, photonic crystal fibers, and random lasers.