Solution (source code)

= Solution

The only positive equilibrium solves $1-y=0$ and $x-1=0$, hence is
$$
\boxed{(x,y)=(1,1)}.
$$
Away from a nullcline,
$$
\frac{dy}{dx}
=\frac{3y(x-1)}{x(1-y)}.
$$
This is separable:
$$
\left(\frac1y-1\right)dy
=3\left(1-\frac1x\right)dx.
$$
Integration gives the <first integral>
$$
\log y-y=3x-3\log x+C.
$$
Thus one convenient conserved quantity is
$$
\boxed{E(x,y)=3x-3\log x+y-\log y}.
$$

Its gradient is
$$
\nabla E=\left(3-\frac3x,\,1-\frac1y\right),
$$
so $(1,1)$ is its only stationary point in the positive quadrant. Its <Hessian matrix> is
$$
\nabla^2E=
\begin{pmatrix}3/x^2&0\\0&1/y^2\end{pmatrix},
$$
which is positive definite everywhere. Hence $(1,1)$ is a strict, indeed global, minimum.

The nearby level curves of $E$ are closed curves surrounding this minimum. Since $E$ is constant along every trajectory, a solution starting on a sufficiently small nearby level set cannot leave the region bounded by a slightly larger level set. This proves <Lyapunov stability> of the equilibrium: solutions initially close to $(1,1)$ remain close for all time.