as big as computationally feasible. Various ways to generate this initial set of wavefunctions for a given $f(x,v)$ are presented in section \[sec:IC\].

*At this stage the precision of the approximation $f(x,v)\to
  \{\psi_n(x)\}$ is controlled by two parameters, $\hbar$ and $N$.*

The wavefunctions $\psi_n(x)$ are now evolved forward in time using the coupled [Schrödinger]{}-Klein-Gordon system of equations $$\begin{array}{rcl}
  i \hbar \partial_t \psi_n &=& - \displaystyle\frac{\hbar^2}{2} \nabla^2_x \psi_n +
  U\psi_n, \Big.\\
  -\displaystyle\frac{1}{c^2}\displaystyle\frac{\partial^2}{\partial t^2}U + \nabla^2_xU &=& 4\pi 
  G\left(\displaystyle\sum_{n=1}^N |\psi_n(x)|^2 - \bar\rho\right). \Big.
\end{array}$$

The integration in time of the [Schrödinger]{}-Klein-Gordon system can be done explicitly using finite differences on a regular grid, as will be described in section \[sec:numerics\].

Controlling your simulation. As the simulation is running you should monitor the following quantities in order to see that the choice of the method does not introduce artefacts. The correction terms $$\label{eq:correction_term}
\sum_{\shortstack{$\scriptstyle r\geq 3$\\$\scriptstyle r~\rm {odd}$}}
\frac{1}{r!}\left(\frac{\hbar}{2i}\right)^{r-1} {\frac{\partial ^r}{\partial x^r}}V {\frac{\partial ^r}{\partial v^r}}f(x,v)$$

should be small when compared to the ones ($v{\frac{\partial f}{\partial x}}$ and ${\frac{\partial U}{\partial x}}{\frac{\partial f}{\partial v}}$) entering the Vlasov equation. Thanks to the $1/r!$ decrease and the smoothness of the gravitational potential $V$ in most cases of interest, computing the first term of this series is generally 