Liouville theorem demonstration

My goal is to create a demonstration of the Liouville theorem in 2D phase space.

I made up an interesting potential energy function $U(x) = (x-4)x^3 + 27$, so that the minimal energy of system is exactly zero.

The Hamiltonian equations are (in dimensionless units):

$$dp/dt = -4(x-3)x^2,\ dx/dt = p$$

The equation for Contour plot is $(1/2)p^2 + U(x) = E$. Now I want to demonstrate that for given initial set of $(x, p)$ in phase plane, forming an area $D(0)$, which is time-evolving, and this area is conserved, i.e, $D(0) = D(t)$. The problem is that this non-linear system of ODEs is not solvable, so I cannot tell for a given initial $(x,p)$, how it will move in phase space after one second, it has to be solved numerically.

The goal is to draw several images -- first with chosen volume $D(0)$ formed by a large number of initial values of (x,p). Next image would be a system evolved in time, say in one second. It'd be nice to see the volume $D(0)$ being distorted due to the Hamiltonian equations and the non-linear force. One approach is 1) to parametrically describe the volume $D(0)$, 2) numerically solve the Hamilton equations for all those points in $D(0)$ to get their position for latter time 3) draw a new image with volume $D(t)$. It should be sufficient to solve the Hamiltonian equations only on the boundary of $D(0)$, which should preserve continuity in time. Another approach is to choose the initial volume $D(0)$, randomly choose, say, a hundred points in this volume, then solve Hamilton equations for all these points until time $t$, and draw these hundred points at their new position.

The problem is, I don't have the slightest idea how to do this practically in Mathematica 8. I'd appreciate any help.

-
I've given you a starting point. If you are looking for a complete solution, people are unlikely to put in the effort to write one, unless you demonstrate that you put in the effort to learn a bit of Mathematica and at least try to solve the problem. – Szabolcs Jul 4 '14 at 17:20
Thank you for the help, I try to use NDsolve for the same equations, but many initial values. – user16320 Jul 4 '14 at 17:35
Be careful not to use E as a symbol for energy in your code. Mathematica reserves this symbol for the base of the natural logarithms. – m_goldberg Jul 4 '14 at 18:46

To solve differential equations numerically, use NDSolve (as a documentation search quickly reveals).
To plot a set of points in 2D, use ListPlot.