A new ansatz for molecular vibronic wave functions based on a superposition of time-dependent Generalized Coherent States is developed and analysed. The methodology is specifically tailored to describe the time evolution of the wave function of a system in which several interacting electronic states are coupled to a bath of harmonic oscillators. The equations of motion for the wave packet parameters are obtained by using the Dirac–Frenkel time-dependent variational principle. The methodology is used to describe the quantum dynamical behavior of a model polaron system and its scaling and convergence properties are discussed and compared with numerically exact results.