We propose a variational approach to the dynamics of the Bose-Hubbard model beyond the mean field approximation. To develop a numerical scheme, we use a discrete overcomplete set of Glauber coherent states and its connection to the generalized coherent states studied in depth by Perelomov [A. Perelomov, Generalized Coherent States and Their Applications, Springer-Verlag (Berlin, 1986)]. The variational equations of motion of the generalized coherent state parameters as well as of the coefficients in an expansion of the wavefunction in terms of those states are derived and solved for many-particle problems with large particle numbers S and increasing mode number M. For M = 6 it is revealed that the number of parameters that have to be propagated is more than one order of magnitude less than in an expansion in terms of Fock states.
