# Verifying that ODE integration generates Theoretical Stationary distribution

I am trying to simulate an ODE, like $\dot{x} = \xi(x)$ that should have a stationary distribution (a la Stat Mech). Assuming that my ODE algorithm generates time samples of my system state $x$ like $x^{(0)}, ..., x^{(t)}, x^{(t+1)}, ..., x^{(n)}$, I would like to verify that the samples generated from the ODE numerically approximately match the theoretical distribution (or at least in metrics like entropy/KL Divergence they match).

I thought about running Monte Carlo and the ODE and trying to compare the generated results, but I don't know how to compare the samples generated as they are generated in a high dimensional phase space. Normally (low dimensions), I would use binning and compute KL divergence.

My motivation is to further convince myself that the results from "On the classical statistical mechanics of non-Hamiltonian systems" by Tuckerman are real.

In general, your type of question would be called a "multivariate goodness of fit test". If $F(x_1,\ldots,x_n)$ is the $n$-dimensional CDF for the theoretical distribution, and the random variables $(X_1,\ldots,X_n)$ are a sample from your ODE, then $$Z_1 = F(X_1), \quad Z_2 = F(X_2\mid X_1), \quad\cdots\quad Z_n = F(X_n\mid X_1,\ldots,X_{n-1})$$ are random variables that are uniformly independently distributed in the unit cube $[0,1]^n$. This is Rosenblatt's transform, similar to the probability integral transform.
So if your system is ergodic and enough time passes between the ODE samples for them to be independent, you can just test that the vectors $z^{(0)},\ldots,z^{(n)}$ are u.i.d. on the unit cube. This in turn can be done, for example, by splitting the cube into bins and using the chi-square test.