  We present numerical methods to solve the Generalized Hartree-Fock theory for
fermionic systems in lattices, both in thermal equilibrium and out of
equilibrium. Specifically, we show how to determine the covariance matrix
corresponding to the Fermionic Gaussian state that optimally approximates the
quantum state of the fermions. The methods apply to relatively large systems,
since their complexity only scales quadratically with the number of lattice
sites. Moreover, they are specially suited to describe inhomogenous systems, as
those typically found in recent experiments with atoms in optical lattices, at
least in the weak interaction regime. As a benchmark, we have applied them to
the two-dimensional Hubbard model on a 10x10 lattice with and without an
external confinement.
