Confined electromagnetic modes strongly couple to collective excitations in ensembles of quantum emitters, producing light-matter hybrid states known as polaritons. Under such conditions, the discrete multilevel spectrum of molecular systems offers an appealing playground for exploring multiphoton processes. This work contrasts predictions from the Tavis-Cummings model in which the material is a collection of two-level systems, with the implications of considering additional energy levels with harmonic and anharmonic structures. We discuss the exact eigenspectrum, up to the second excitation manifold, of an arbitrary number N of oscillators collectively coupled to a single cavity mode in the rotating-wave approximation. Elaborating on our group-theoretic approach [New J. Phys. 23, 063081 (2021)], we simplify the brute-force diagonalization of N2 × N2 Hamiltonians to the eigendecomposition of, at most, 4 × 4 matrices for arbitrary N. We thoroughly discuss the eigenstates and the consequences of weak and strong anharmonicities. Furthermore, we find resonant conditions between bipolaritons and anharmonic transitions where two-photon absorption can be enhanced. Finally, we conclude that energy shifts in the polaritonic states induced by anharmonicities become negligible for large N. Thus, calculations with a single or few quantum emitters qualitatively fail to represent the nonlinear optical response of the collective strong coupling regime. Our work highlights the rich physics of multilevel anharmonic systems coupled to cavities absent in standard models of quantum optics. We also provide concise tabulated expressions for eigenfrequencies and transition amplitudes, which should serve as a reference for future spectroscopic studies of molecular polaritons.