I am tasked to find the Hamiltonian of a system which comprises of a body of mass $m$ that moves on the surface defined by $z = x(x-1) + y(y-1)$. I have computed the Lagrangian to this problem and found, $L = m\dot{x}^2(2x^2 -2x +1)+m\dot{y}^2(2y^2 -2y +1)+2m\dot{x}\dot{y}(4xy - 2x -2y +1) + mg(x(x-1) + y(y-1)).$
I have computed the $p_x$ and $p_y$ and found them to be, $p_x = 2m\dot{x}(2x^2 - 2x +1) + 2m\dot{y}(4xy - 2x -2y+1)$ and $p_y = 2m\dot{y}(2y^2 - 2y +1) + 2m\dot{x}(4xy - 2x -2y+1)$; $p_x = \dfrac{\partial L}{\partial \dot{x}}$ and $p_y = \dfrac{\partial L}{\partial \dot{y}}$.
I can compute the values of $\dot{x}$ and $\dot{y}$ from here and then plugging those values in $\dot{x}p_x + \dot{y}p_y - L \tag{$\star$},$ I will obtain the required Hamiltonian.
The problem is,
$\dot{x} = \dfrac{p_y(4xy - 2x - 2y +1) - (2y^2 -2y +1)p_x}{2m((12(x^2y^2 - x^2y - xy^2 + xy) + 2(x^2 + y^2 -x-y))}$. Which is a huge expression and plugging that into $(\star)$ is getting very cumbersome.
So, I want to know if there is any other way to obtain the Hamiltonian or any manipulation that I can do to make my task easier.
Thanks!