  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.
