I am looking for a 2D ODE of the form:
$\begin{array}{rl} \dot x&= f(x,y) \\ \dot y &= g(x,y) \end{array}$
that satisfies the following:
1- there exists a function $H(x,y)$ such that $f(x,y)=-\frac{\partial H}{\partial x}, g(x,y)=-\frac{\partial H}{\partial y}$.
2- The ODE is positive, which means the positive quadrant is forward invariant. This means that if $x(0)>0,y(0)>0$ then the ODE solution satisfies $x(t)>0,y(t)>0$ for all t. This is usually verified if: $\forall y>0, f(0,y)>0$ and $\forall x>0, g(x,0)>0$
3- there are at least two steady states of the ODE which are positive. ( a steady state $(x_e,y_e)$ is positive if $x_e>0,y_e>0$).
4- Each of the aforementioned steady states is stable, i.e., the eigenvalues of the Jacobian have negative real parts when evaluated at each of them.
The function $f(u,v) = e^{-4v} - 2 u^4 - 4 u^2 e^{-v}$ has local minima at $(u,v)=(\pm 1,0)$ and no other stationary points, so we can shift and rotate the coordinate system to obtain a function of the kind you seek, with local minima at $(x,y)=(3,1)$ and $(1,3)$: $$ H(x,y) = f \bigl( \tfrac{x-y}{2}, \tfrac{x+y}{2}-2 \bigr) = e^{8-2x-2y} - \tfrac18 (x-y)^2 - (x-y)^2 e^{2-\tfrac12(x+y)} . $$ Here's the resulting phase portrait for $(\dot x,\dot y) = - \nabla H$ in the positive quadrant, with the nullclines in red and orange: