A reliable and efficient computation of the entire single-particle spectrum of correlated molecules is an outstanding challenge in the field of quantum chemistry, with standard density functional theory approaches often giving an inadequate description of excitation energies and gaps. In this work, we expand upon a recently-introduced approach which relies on a fully self-consistent many-body perturbation theory, coupled to a non-perturbative truncation of the effective dynamics at each step. We show that this yields a low-scaling and accurate method across a diverse benchmark test set, capable of treating moderate levels of strong correlation effects, and detail an efficient implementation for applications up to $sim1000$ orbitals on parallel resources. We then use this method to characterise the spectral properties of the artemisinin anti-malarial drug molecule, resolving discrepancies in previous works concerning the active sites of the lowest energy fundamental excitations of the system.