We present a separable approach to linear response in metal clusters, employing a sum of separable terms of local operators. The coefficients of the separable force are unambiguously derived from the full residual interaction of the system. Comparing with results from the full RPA treatment, we find that the separable approach provides a reliable and detailed description of the strength function at much less expense. As a first application, we investigate the strength distribution for the surface plasmon resonance in neutral and singly ionized, spherical Na clusters over a wide range of sizes, up to N = 5108. There is surprisingly much structure in the resonance spectra even for the very large clusters. The asymptotic trends for peak position and FWHM of neutral clusters in the steep jellium model agree perfectly well with previously derived analytical results. Taking into account a finite surface in the jellium model and using charged clusters, the asymptotic slopes are substantially increased. In both cases, the asymptotic trends develop only very late, i.e. for cluster sizes safely above N = 1000.