To a large extent, our knowledge of the photoinduced dynamics of molecular systems at the atomic level is shaped by nonlinear femtosecond spectroscopy. Traditionally, nonlinear femtosecond spectroscopy is an ensemble spectroscopy, performed on ensembles of identical chromophores in the gas phase or in the liquid phase. Modern femtosecond spectroscopy comprises a set of various third-order or higher-order technique including, for example, fluorescence up conversion, transient absorption, and photon echo spectroscopy. Very recently, the portfolio of femtosecond techniques has been extended towards single-molecule detection by the development of fluorescence-detected double-pump single-molecule spectroscopy. With this technique, a temporal resolution of about ten femtoseconds has been achieved. The technique permits the real-time monitoring of not only electronic populations, but also of electronic and vibrational coherences for individual molecules. Given the complexity of systems which are of interest in current applications, the simulation of femtosecond signals is a major challenge for theory. We have developed general theories for the simulation of ensemble as well as single molecule spectroscopy beyond the weak-field limit. The polarization obtained by this theory reduces to the nonlinear response function formalism in the limit of weak field-matter interaction. For intermediate and strong pulse strengths, this theory takes into account all relevant higher-order contributions beyond the weak-field limit.  We have applied our theory to simulate ensemble as well as single molecule spectroscopy of various light harvesting systems, revealing intricate energy transfer processes in these systems.