Phase-space representations are a family of methods for dynamics of both bosonic , fermionic systems, that work by mapping the system's density matrix to a quasiprobability density and the Liouville-von Neumann equation of the Hamiltonian to a corresponding density differential equation for the probability. We investigate here the accuracy and the computational efficiency of one approximate phase-space representation, called the fermionic truncated Wigner approximation (fTWA), applied to the Fermi-Hubbard model. On a many-body 2D system, with hopping strength and Coulomb U tuned to represent the electronic structure of graphene, the method is found to be able to capture the time evolution of first-order (site occupation) and second-order (correlation functions) moments significantly better than the mean-field, Hartree-Fock method. The fTWA was also compared to results from the exact diagonalization method for smaller systems , in general the agreement was found to be good. The fully parallel computational requirement of fTWA scales in the same order as the Hartree-Fock method, and the largest system considered here contained 198 lattice sites.