We adapted a recently proposed time-domain framework to characterize the optical response of interacting electronic systems in order to expedite its computation without compromise in quantitative or qualitative accuracy at the microscopic level. With reliable parameterizations of Hamiltonians and interactions, our formulation allows for increased economy and flexibility in calculating the optical response functions to fields of arbitrary temporal shape and strength. For example, the computation of high-harmonic susceptibilities to arbitrary order becomes straightforward within a unified scheme that natively takes into account excitonic effects, as well as deviations of the electronic system from equilibrium under a strong field. Given that two-dimensional semiconductors are currently of much interest for their strong optical nonlinearities, largely defined by excitons, we demonstrate the approach by computing the frequency-dependent susceptibilities of monolayer MoS2 and hexagonal boron nitride up to the third-harmonic. In the latter, a two-band model brings further insight on the role of intra-band transitions and the nonequilibrium state of the system when computing even-order response, like the second-harmonic susceptibility. Being grounded on a generic non-equilibrium many-body perturbation theory, this framework allows extensions to handle more generic interaction models or the realistic description of electronic processes taking place at ultrafast time scales.