I am trying to show that the function $$Q(x,y)=be(x-\bar{x})^2+2ce(x-\bar{x})(y-\bar{y}) +cf(y-\bar{y})^2$$ is a Lyapunov function for the competition equations: $$\dot{x}=x(a-bx-cy)$$ $$\dot{y}=y(d-ex-fy).$$ I can find: $$(\bar{x},\bar{y})=\Big(\frac{af-cd}{ce-bf},\frac{bd-ae}{bf-ce}\Big). $$
I know that for a function, $V(x,y)$ to be Lyapunov, it must be $$V(x,y) \geq 0, \; \forall x,y \in \mathbb{R},$$ how can this be seen? and that $$\dot{V} \leq 0 \; \forall x,y \in \mathbb{R}.$$
Am I missing a trick? I am just getting a lot of messy algebra. I cannot seem to get $Q$ into any kind of addition of squared quantities, nor can i figure out a way to make the derivative negative for all point in the plane. $$Q(x,y)=be(x-\bar{x})^2+2ce(x-\bar{x})(y-\bar{y}) +cf(y-\bar{y})^2$$ $$\dot{Q}(x,y)=2\Big[\big(be(x-\bar{x})+ce(y-\bar{y})\big)\dot{x}+\big(ce(x-\bar{x}) +cf(y-\bar{y})\big)\dot{y}\Big]$$
The dynamics can also be written as
$$ \dot{z} = \text{diag}(z) (v - M\,z), \tag{1} $$
with
$$ z = \begin{bmatrix} x \\ y \end{bmatrix}, \quad v = \begin{bmatrix} a \\ d \end{bmatrix}, \quad M = \begin{bmatrix} b & c \\ e & f \end{bmatrix}. $$
From this it is also easier to see what the equilibria ($\dot{z}=0$) are, namely $\bar{z}_1=0$ or $\bar{z}_2 = M^{-1}v$. The second equilibrium actually gives
\begin{align} \bar{x}_2 &= \frac{a\,f - c\,d}{b\,f - c\,e}, \tag{2a} \\ \bar{y}_2 &= \frac{b\,d - a\,e}{b\,f - c\,e}, \tag{2b} \end{align}
which differs in a minus sign in front of $\bar{x}$ compared with your equilibrium. The candidate Lyapunov function can also be written into a "more compact" form
\begin{align} Q(x,y) &= \begin{bmatrix} x - \bar{x} \\ y - \bar{y} \end{bmatrix}^\top \begin{bmatrix} be & ce \\ ce & cf \end{bmatrix} \begin{bmatrix} x - \bar{x} \\ y - \bar{y} \end{bmatrix} \tag{3a} \\ &= (z - \bar{z})^\top P\,(z - \bar{z}) \tag{3b} \end{align}
with
$$ P = \begin{bmatrix} b\,e & c\,e \\ c\,e & c\,f \end{bmatrix} \tag{4}. $$
By defining error coordinates as $\varepsilon = z - \bar{z}_2$, then its dynamics can be written as
\begin{align} \dot{\varepsilon} &= \text{diag}(z) (v - M\,z) \tag{5a} \\ &= \text{diag}(z) (v - M\,(\varepsilon + M^{-1}v)) \tag{5b} \\ &= -\text{diag}(z)\,M\,\varepsilon \tag{5c} \end{align}
This allows one to simplify the candidate Lyapunov function even further to
$$ Q(x,y) = \varepsilon^\top P\,\varepsilon \tag{6}. $$
The constraint that this should satisfy $Q(x,y) \geq 0\ \forall\,x,y\in\mathbb{R}$ is equivalent to that all the eigenvalues of $P$ need to be non-negative. However, it can be noted that $Q(x,y)=0$ would also satisfy this requirement. So instead I will use that $Q(x,y) > 0\ \forall\,\varepsilon\neq0$. This is then equivalent to that all the eigenvalues of $P$ need to be positive. The expression for the eigenvalues are a bit cumbersome. Since $P\in\mathbb{R}^{2\times 2}$ we have $\det(P) = \lambda_1\,\lambda_2$ and $\text{tr}(P) = \lambda_1 + \lambda_2$, with $\lambda_1$ and $\lambda_2$ the two eigenvalues of $P$. So instead of looking at the eigenvalues directly we can also use that $\det(P) > 0$ and $\text{tr}(P) > 0$
\begin{align} \det(P) &= c\,e\,(b\,f - c\,e) > 0, \tag{7a} \\ \text{tr}(P) &= b\,e + c\,f > 0. \tag{7b} \end{align}
By using $(5c)$ and $(6)$ the time derivative of the candidate Lyapunov function can be expressed as
\begin{align} \dot{Q}(x,y) &= \dot{\varepsilon}^\top P\,\varepsilon + \varepsilon^\top P\,\dot{\varepsilon} \tag{8a} \\ &= -\left(\text{diag}(z)\,M\,\varepsilon\right)^\top P\,\varepsilon - \varepsilon^\top P\,\text{diag}(z)\,M\,\varepsilon \tag{8b} \\ &= -\varepsilon^\top \left(M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M\right) \varepsilon \tag{8c} \end{align}
Checking whether this satisfies the constraint that $\dot{Q}(x,y) \leq 0\ \forall\,x,y\in\mathbb{R}$ becomes a bit more cumbersome. Namely one could look at whether the eigenvalues of $M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M$ are always non-negative, but this would be a function of $z$ and thus also $\varepsilon$ so even if an eigenvalue is negative for a given $x$ and $y$ it might still be that $\dot{Q}(x,y) \leq 0$. Only if both eigenvalues are negative then $\dot{Q}(x,y) \leq 0$ would be false. So stability is unsure if only the determinant or the trace (so not both at the same time) is positive. If they both are non-negative it is stable, but when evaluating the determinant and trace it can be seen that that is not the case
\begin{align} \det(M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M) &= 4\,c\,e\,x\,y\,(b\,f - c\,e)^2 > 0, \tag{9a} \\ \text{tr}(M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M) &= 2 \left(e\,x\,(b^2 + c^2) + c\,y\,(e^2 + f^2)\right) > 0, \tag{9b} \end{align}
since for example depending on the sign of $c\,e$ a different pair of opposing quadrants for $x$ $y$ will yield that the determinant is negative. But maybe since you are dealing with a competitive Lotka–Volterra equation you are only interested in $x,y\in\mathbb{R}_{\geq0}$. But I am going to leave it here, so hopefully this answer has given you enough insights and tricks to help to proceed.