We present a variational solution of the time-dependent Schrodinger equation formed from a restricted superposition of frozen Gaussian wavepackets. The trial function is comprised of a set of frozen Gaussian wavepackets and a set of three time-dependent variational parameters. The trial wavefunction is subjected to the McLachlan variational principle which leads to a set of equations for the optimal time-evolution of the variational parameters. We present numerical results for the time-evolution of the trial wavefunction on a single-dimensional anharmonic potential energy and compare to the results to the exact time-evolution. Finally, we derive a multi-dimensional generalization of the algorithm. (c) 2005 Elsevier B.V. All rights reserved.