  We present a straightforward implementation scheme for solving the time
dependent Schr\"odinger equation for systems described by the Hubbard
Hamiltonian with time dependent hoppings. The computations can be performed for
clusters of up to 14 sites with in principle general geometry. For the time
evolution, we use the exponential midpoint rule, where the exponentials are
computed via a Krylov subspace method, which only uses matrix-vector
multiplication. The presented implementation uses standard libraries for
constructing sparse matrices and for linear algebra therefore the approach is
easy to use on both desktop computer and computational cluster. We apply the
method to calculate time evolution of double occupation and nonequilibrium
spectral function of a photo-excited Mott-insulator. The results show that not
only the double occupancy increases due to creation of electron-hole pairs but
also the Mott gap becomes partially filled.
