$$\dot{x} = x - x^2 - xy$$
$$\dot{y} = y - y^2 - xy$$
Is there a general method to find a Lyapunov function? Why do I feel like finding a Lyapunov function is like shooting in the dark and luck-dependent? I tried to guess some function with $x$ and $y$, but it doesn't work. I also tried Mathematica, and apparently, this is too complicated for Mathematica to solve.
A side-note question: if the general method to find a Lyapunov function is currently not available, is that because no one is smart enough to come up with the general method? Or is it due to its nature that it's "impossible" to have the general method? And has anyone proved its impossibility?
Note: it's not a duplicate question because only my side-note question is a duplicate of the other question. My first question is tailored to this specific system, which is unique.
First an observation on the direct solution of the system. Add both equations to get $$ \frac{d}{dt}(x+y)=(x+y)-(x+y)^2 $$ which is a logistic equation for $u=x+y$ with a stable point at $x+y=1$ and an unstable at $x+y=0$, $$ \dot u=u-u^2\text{ or }\frac{du^{-1}}{dt}=1-u^{-1}\implies 1-u^{-1}=(1-u_0^{-1})e^{-t}. $$ Then building on that solution consider the difference of the original equations $$ \frac{d}{dt}(x-y)=(x-y)-(x^2-y^2)=(x-y)(1-x-y) $$ so that $v=x-y$ satisfies $$ \frac{\dot v}{v}=1-u=\frac{\dot u}{u}\implies v(t)=\frac{v_0}{u_0}u(t) $$ Any solution starting in $(x_0,y_0)$ with $x_0+y_0>0$ converges to $$(x_*,y_*)=\frac{(x_0,y_0)}{x_0+y_0}$$ along a straight line, all solutions starting from the other half-plane $x_0+y_0\le 0$ diverge exponentially. All that along straight lines originating at the origin, which means that $(0,0)$ is an unstable point.
The "standard" Lyapunov function $V=x^2+y^2$ confirms that, as $$ \dot V=2V(1-x-y)\ge 2V(1-\sqrt{2V}) $$ which is positive for $\|(x,y)\|^2=V\le\frac12$ so that the origin is repelling.