By employing the Dirac-Frenkel time-dependent variational principle, we study the dynamical properties of the Holstein molecular crystal model with diagonal and off diagonal exciton-phonon coupling. A linear combination of the Davydov D-1 (D-2) ansatz, referred to as the "multi-D-1, ansatz" ("multi-D-2 ansatz"), is used as the trial state with enhanced accuracy but without sacrificing efficiency. The time evolution of the exciton probability is found to be in perfect agreement with that of the hierarchy equations of motion, demonstrating the promise the multiple Davydov trial states hold as an efficient, robust description of dynamics of complex quantum systems. In addition to the linear absorption spectra computed for both diagonal and off-diagonal cases, for the first time, 2D spectra have been calculated for systems with off-diagonal exciton-phonon coupling by employing the multiple D-2 ansatz to compute the nonlinear response function, testifying to the great potential of the multiple D-2 ansatz for fast, accurate implementation of multidimensional spectroscopy. It is found that the signal exhibits a single peak for weak off-diagonal coupling, while a vibronic multipeak structure appears for strong off-diagonal coupling.