Vibrationally resolved electronic spectroscopy provides important information on the structure and dynamics of polyatomic molecules. On-the-fly ab initio implementation [1,2,3] of the thawed Gaussian approximation (TGA) [4] has been successful in describing such spectra within the Condon approximation but cannot capture non-Condon effects, tunneling or wavepacket splitting. While the extended TGA [5,6] captures first-order non-Condon effects, Hagedorn’s wavepackets [7,8], i.e. superposition of states, each of which is the Gaussian multiplied by a cleverly chosen polynomial, enable the propagation of non-Condon wavepackets of arbitrary shapes. I will demonstrate that a variational implementation [9] of the Hagedorn wavepackets can capture tunneling and wavepacket splitting, as well as their effects on spectra. To demonstrate that the method is not restricted to small model systems, I will combine our recently developed algorithm for an algebraic evaluation of the overlap of Hagedorn wavepackets with on-the-fly ab initio semiclassical dynamics to compute single vibronic level fluorescence spectra of polyatomic molecules.