We present a real-time propagation method for computing linear and nonlinear optical properties of molecules based on the Bethe–Salpeter equation. The method follows the time evolution of the one-particle density matrix under an external electric field. We include electron–electron interaction effects through a self-energy model based on the screened exchange approximation. Quasiparticle energies are taken from a prior GW calculation to construct the effective single-particle Hamiltonian, and we represent all operators and wave functions in an atom-centered Gaussian basis. We benchmark the accuracy of the real-time propagation against the standard linear-response Bethe–Salpeter equation by using a set of organic molecules. We find very good agreement when computing linear-response isotropic polarizability spectra from both approaches with a mean absolute deviation of 30 meV in peak positions. Beyond linear response, we simulate second harmonic generation and optical rectification in a noncentrosymmetric molecule. We foresee broad applicability of real-time propagation based on the Bethe–Salpeter equation for the study of linear and nonlinear optical properties of molecules, as the method has a computational cost similar to that of time-dependent density functional theory with hybrid functionals.