You can add the bouncing behavior using a mirroring method, where we relax the system and instead of bouncing the ball one changes the sign of the potential. It can be written as
$H=\frac{p^2}{2m}+mg|y|$
where $|y|=y$ if $y>=0$ and $|y|=-y$ if $y<0$. Potential then looks like (for $m,g=1$)
Equation of motion becomes
$\ddot y(t)=-mg\,sgn(y)$,
where $sgn(y)$ is the sign function of y (a Heaviside-theta-like function). To allow equlilibrium at $y=0$ to exist, we can define $sgn(0)\equiv0$.
Writing this in Wolfram Alpha (here, with $m=-g=1$ for simplicity and y=q) gives the following graphs related to explicit solution $y(t)$ and phase diagram trajectories:
Inserting the Hamiltonian explicitly gives a better look at phase space appearance
Explicit solution to $y(t)$ has an iterative character but in essence is similar to the free-fall solution (it only contains an extra $sgn$ function in the $t^2$ term):
$y(t)=y_0+v_0t-\frac g2sgn(y)t^2$.
And, finally, to recover the actual bouncing ball behavior just take the absolute value of $y(t)$,
$|y(t)|$.
This model is very interesting in that it takes a very simple and linear (in fact constant) model and turns it into a highly non-linear system (iterative, as shown), by adding an impulse at $y=0$ (the bouncing behavior).