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.
You can solve for the trajectories directly. Observe,
\begin{align}
\frac{dx}{dy} &= \frac{\dot{x}}{\dot{y}}\\
&= \frac{x + y\cos(y)}{-y}
\end{align}
By Mathematica, this has the general solution
\begin{align}
x &= \frac{C - \cos(y)}{y} - \sin(y) \,,\quad C \in \mathbb{R}
\end{align}
We need $\lim_{y \to 0} x(y) = 0$, so we infer that $C = 1$. Thus, we have the equation for the unstable manifold:
\begin{align}
\boxed{x = \frac{1 - \cos(y)}{y} - \sin(y) }
\end{align}
To verify that this is the case, we can view a plot of the phase portrait:
If you want to play around with this, here's some Mathematica code:
sol = DSolve[{x'[y] == (x[y] + y*Cos[y])/-y, x[0] == 0}, x[y],
y][[1]] // FullSimplify
xmin = -4; xmax = 4; ymin = -3; ymax = 3;
cplot = ContourPlot[
x == (x[y] /. sol), {x, xmin, xmax}, {y, ymin, ymax},
ContourStyle -> {Black}];
splot = StreamPlot[{x + y*Cos[y], -y}, {x, xmin, xmax}, {y, ymin,
ymax}, StreamColorFunction -> "Rainbow"];
Manipulate[
Show[splot, cplot,
ParametricPlot[
Evaluate[
First[{x[t], y[t]} /.
NDSolve[{x'[t] == x[t] + y[t]*Cos[y[t]], y'[t] == -y[t],
Thread[{x[0], y[0]} == point]}, {x, y}, {t, 0, T}]]], {t, 0,
T}, PlotStyle -> Red]], {{T, 20}, 1, 100}, {{point, {3, 0}},
Locator}, SaveDefinitions -> True]
Best Answer
We have $$ \frac12\frac{d}{dt}(x^2+y^2)=-x^2 + x y - y^2/4 - y^4 $$ and the matrix in the quadratic form $$ -x^2 + x y - y^2/4=\begin{pmatrix}x& y\end{pmatrix}\begin{pmatrix}-1 & 1/2 \\ 1/2 &-1/4\end{pmatrix} \begin{pmatrix} x\\ y\end{pmatrix} $$ has eigenvalues $-5/4$ and $0$, the latter with eigenvector $(1/2,1)$. This means that $$-x^2 + x y - y^2/4\le0$$ with equality if and only if $y=2x$. Therefore, $$ \frac12\frac{d}{dt}(x^2+y^2)< - y^4 $$ for $y\ne 2x$, and $$ \frac12\frac{d}{dt}(x^2+y^2)\le - y^4 $$ for $y= 2x$. Summing up we have $$ \frac12\frac{d}{dt}(x^2+y^2)<0 $$ outside the origin.