Consider some dynamical system $\dot{\textbf{X}}=F(\textbf{X})$ where $\textbf{X}$ is discretized along a 1-dimensional spatial coordinate $\textbf{x}=(x_1,\dots,x_N)^T$. Let $\rho(\textbf{X},t)$ be the probability density function, then the Liouville equation describes the evolution of $\rho(\textbf{X},t)$ given by,
$\begin{align} \frac{\partial \rho}{\partial t}=-\sum_{j=0}^{N}\frac{\partial}{\partial X_j}(F(\textbf{X})\rho(\textbf{X},t)) \end{align}$
Here is where I may be confusing myself. If we want to solve this numerically then we need to discretize the phase space. I will assume that $X_j \in \mathbb{R}^L$ for each $j$. This means that $\rho$ is a lattice in an ultra-high dimensional phase space of size $L^N$. The Liouville equation is therefore out of reach to solve numerical for even modestly sized spatial grids.
My question is, am I thinking about this correctly? I spoke to someone who said that I am double discretizing, which is a mistake. However, this is the way that I see the problem so I am looking for either a verification of this method or for someone to point out how I am making an error.