Finite lattice models are a prototype for strongly correlated quantum systems
and capture essential properties of condensed matter systems. With the dramatic
progress in ultracold atoms in optical lattices, finite fermionic Hubbard
systems have become directly accessible in experiments, including their
ultrafast dynamics far from equilibrium. Here, we present a theoretical
approach that is able to treat these dynamics in any dimension and fully
includes inhomogeneity effects. The method consists in stochastic sampling of
mean-field trajectories and is found to be more accurate and efficient than
current nonequilibrium Green functions approaches. This is demonstrated for
Hubbard clusters with up to 512 particles in one, two and three dimensions