Previously, the non-thermal emission from galaxies has only been modeled with single-zone models which is insufficient to explain a multitude of new, spatially resolved multi-messenger data of cosmic ray (CR) spectra, at gamma-rays and in the radio. Instead, we perform high-resolution magneto-hydrodynamic simulations of galaxies using the moving mesh code AREPO with selfconsistent CR physics. In post-processing, we calculate steady-state spectra of CRs including all relevant cooling and escape losses. Consistent with Voyager-1 and AMS-02 data, our models show a turn-over of proton spectra below GeV energies due to Coulomb interactions so that electrons start to dominate the total particle spectra and match the shape of the positron fraction up to 10 GeV. Furthermore, from our CR spectra, we calculate multi-frequency spectra, from the radio up to the TeV energy regime, due to all non-thermal emission processes, i.e. synchrotron, bremsstrahlung, inverse Compton (IC) emission and gamma-ray emission from neutral pion decay. This allows us to produce detailed emission maps, luminosities and spectra of our simulated galaxies, that range from dwarfs to Milk-Way analogues to starburst galaxies, at different evolutionary stages. Within our simulations, we can successfully reproduce the observed far infrared (FIR)-gamma-ray and FIR-radio relations. We find that highly SF galaxies are close to the calorimetric limit and hence, their gamma-ray emission is dominated by neutral pion decay. However, in low SF galaxies, escape losses due to diffusion steepen the spectra and in turn, an increasing contribution from IC emission is needed to reproduce the observed gamma-ray spectra.