  We introduce a new second order in time Besse-type relaxation scheme for
approximating solutions of the Schr\"odinger-Poisson system. More specifically,
we use the Crank-Nicolson scheme as a time stepping mechanism, the standard
conforming finite element method for the spatial discretization whilst the
nonlinearity is handled by means of a relaxation approach similar to the one
introduced by Besse for the nonlinear Schr\"odinger equation \cite{Besse}. We
prove that discrete versions of the system's conservation laws hold and we
conclude by presenting some numerical experiments, including an example from
cosmology, that demonstrate the effectiveness and robustness of the new scheme.
