The modeling of vibrations in optical spectra relies heavily on the simplifications brought about by using harmonic oscillators. However, realistic molecular systems can deviate substantially from this description. We develop two methods which show that the extension to arbitrarily shaped potential energy surfaces is not only straightforward, but also efficient. These methods are applied to an electronic two-level system with potential energy surfaces of polynomial form and used to study anharmonic features such as the zero-phonon line shape and mirror-symmetry breaking between absorption and fluorescence spectra. The first method, which constructs vibrational wave functions as linear combinations of the harmonic oscillator wave functions, is shown to be extremely robust and can handle large anharmonicities. The second method uses the cumulant expansion, which is readily solved, even at high orders, thanks to an ideally suited matrix theorem.