systems with underlying bosonic structure, exploiting both quantum interference and an intrinsically large Hilbert space. Here, we experimentally realize a bosonic superconducting processor that combines arbitrary state preparation, a complete set of Gaussian operations, plus an essential non-Gaussian resource – a novel single-shot photon number resolving measurement scheme – all in one device. We utilize these controls to simulate the bosonic problem of molecular vibronic spectra, extracting the corresponding Franck-Condon factors for photoelectron processes in H2O, O3, NO2, and SO2. Our results demonstrate the versatile capabilities of the circuit QED platform, which can be extended to include non-Gaussian operations for simulating an even wider class of bosonic systems.