Existence and uniqueness of homoclinic orbit

Let me explain the idea from my previous answer with few more additional details.

As it was mentioned in the question, the system has a smooth first integral $H(\mathbf{x})$ — a (locally non-constant) function that is constant along the system's trajectories. The key observation is that if such system has a trajectory $\gamma(t)$ homoclinic to a saddle $p$, then $H(p) \equiv H(\gamma(t))$. To prove that we can pick a sequence of moments $t_i \rightarrow +\infty$. By continuity, since $\gamma(t_i) \rightarrow p$ and $H(\mathbf{x})$ is smooth, then $H(\gamma(t)) \rightarrow H(p)$. But since $H(\mathbf{x})$ is constant at any point of $\gamma(t)$, we have that $H(\gamma(t)) \equiv H(p)$.

This gives us the following method to find all homoclinic trajectories for a 2D system of differential equations. First, find all saddle equilibria. Second, find the level sets of $H(\mathbf{x})$ that contain these saddles. What we have proven before means that a homoclinic trajectory must lie in the same level set as the saddle, to which it is homoclinic. Just take a look at level sets after that and you can find homoclinic (and even heteroclinic) trajectories this way.