Solution (source code)

= Solution

Since $p_y=0$ initially and remains constant, the motion stays in a plane of constant $y$. Setting $P=\sqrt{2mH_0}$ and $p=p_{z0}$ gives
$$
\boxed{\dot z=\frac{p+ax^2}{m},\qquad
\dot x=\pm\frac{\sqrt{P^2-(p+ax^2)^2}}m.}
$$
The allowed interval obeys $|p+ax^2|\le P$. The sign of $\dot x$ is set by the initial direction and changes at a turning point. The specified $x_0$ fixes the starting <position> on the <energy> orbit, but the data do not specify which sign of the initial transverse <velocity> is chosen.