  We introduce a new structure preserving, second order in time relaxation-type
scheme for approximating solutions of the Schr\"odinger-Poisson system. More
specifically, we use the Crank-Nicolson scheme as a time stepping mechanism,
whilst the nonlinearity is handled by means of a relaxation approach in the
spirit of \cite{Besse, KK} for the nonlinear Schr\"odinger equation. For the
spatial discretisation we use the standard conforming finite element scheme.
The resulting scheme is explicit with respect to the nonlinearity, i.e. it
requires the solution of a linear system for each time-step, and satisfies
discrete versions of the system's mass conservation and energy balance laws for
constant meshes. The scheme is seen to be second order in time. We conclude by
presenting some numerical experiments, including an example from cosmology and
an example with variable time-steps which demonstrate the effectiveness and
robustness of the new scheme.
