The quantum mechanical description of manyelectron dynamics in molecules driven by short laser pulses is at the heart of theoretical attochemistry. In addition to the formidable time-dependent electronic structure problem, the field faces the challenge that nuclear motion, ideally also treated quantum mechanically, may not be negligible, but scales enormously in effort. As a consequence, most first-principles calculations on ultrafast electron dynamics in molecules are done within the fixed-nuclear approximation. For laser-pulse excitation in H2+, where an "exact" treatment of the coupled nuclear-electron dynamics is possible, it has been shown that nuclear motion can have a nonnegligible impact on high harmonic generation (HHG) spectra (Witzorky et al., J. Chem. Theor. Comput. 2021, 17, 73537365). It is not so clear, however, how to include (quantum) nuclear motion also for more complicated molecules, with more electrons and/or nuclei, in particular when the electronic structure is described by correlated, multistate wavefunction methods such as the time-dependent configuration interaction (TD-CI). In this work, we suggest a scheme in which the Born-Oppenheimer potential energy surfaces of a molecule are approximated by model potentials (harmonic and asymptotic, as an expansion in 1/R), obtained from only a few ab initio calculations, with the prospect to treat complex molecular systems. The method is tested successfully for HHG by few-cycle laser pulses for the "exact" H2+ reference. It is then applied for diatomic molecules with more electrons and for a two-dimensional model of the water molecule using TD-CIS (S = single) for the electronic structure part.