How to find a Lyapunov function? $$\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. 
 A: 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.
A: Actually there is no straightforward method to calculate Lyapunov function but you should care that we compute Lyapunov function to show the stability of a dynamic system about its equilibrium point, where the derivative of the variables are equal to zero. In this case  $(x,y)=(1,0) , (0,0)$. You should try to find a Lapunov function
A: Excuse me. You should try to find a Lyapunov function if the dynamic system is stable around one of these equilibrium points. By linearixation the system around these two equilibrium points we have these two coefficient matrixes:
$$
\begin{matrix}
1 & 0 \\
0 & 1 \\
\end{matrix}
$$
$$
\begin{matrix}
1 & 0 \\
0 & 1 \\
\end{matrix}
$$
If you calculate the eigenvalues of these matrixes simply you will find out that both of them have an eigenvalue in the right hand side of the complex numbers surface. So both of the equilibrium points are unstable in the view of the linearized system. Since the linearized system is unstable, the primary system will be unstable too. So the system is unstable in its both equilibrium point and we shoul not try to prove its stability in its equilibriums. Indeed you shoul try to prove its unstability by using method such that presented here or a method such as Cheatev theorem. 
A: I can not add any comment but this is not true to decide about the stability of each part of a summation by knowing that the sum is stable.Suppose that x=f(t) tend to infinity so y=1-f(t) tend to infinity too. but x+y=1 is stable.In addition when we want to claim about stability we should not consider any constraint over the initial value. In the view of the control theory both of the equilibrium point of this system are unstable.
