Let me repeat, with suitable modifications, my answer to the OP's earlier question.
Assume that all solutions starting in $\mathbb{R}^2_{+} := \{\, (x, y) : x \ge 0,\ y \ge 0 \,\}$ are bounded for $t > 0$ (this would require a separate proof). It follows then that for any such solution, its domain contains $[0,\infty)$ and its $\omega$-limit set is compact and nonempty.
Let $L$ stand for the $ω$-limit set of some point, $(x_0,y_0)$, sufficiently close to the unstable focus $(1,3)$. By the Poincaré–Bendixson theorem, as there are finitely many equilibria, $L$ is either a limit cycle, or a heteroclinic cycle, (EDIT: or a homoclinic loop), or an equilibrium.
There are two equilibria, $(4,0)$ and $(1,3)$. The first of them is an unstable focus, so it cannot belong to any heteroclinic cycle (because a heteroclinic cycle (EDIT: or a homoclinic loop) must contain an equilibrium that is an $\omega$-limit point for some other point). Consequently, there are no heteroclinic cycles (EDIT: or homoclinic loops) at all.
So, $L$ is either a periodic orbit, or equals $\{(4,0)\}$. We proceed now to excluding the latter. The linearization of the vector field at $(4,0)$ has matrix
$$
\begin{bmatrix}
-1 & -\frac25
\\
0 & 3
\end{bmatrix}.
$$
There are two real eigenvalues, $-1$ and $3$, of opposite signs, so $(4,0)$ is a hyperbolic saddle. Its stable manifold is tangent to an eigenvector corresponding to $-1$, that is, to $(1,0)$. I claim that the stable manifold
is the $x$-axis, minus $(4,0)$. Indeed, on the $x$-axis we have $\dot{x} = 4 - x$, $\dot{y} \equiv 0$, so for any $(x_1,0)$ we have $\omega((x_1,0)) = \{(4,0)\}$. Now, for a hyperbolic saddle its stable manifold is just the set of those points whose (unique) $\omega$-limit point is the saddle. Hence, if $L = \{(4,0)\}$ then the positive semitrajectory of $(x_0, y_0)$ must belong to the $x$-axis, which contradicts the uniqueness of the initial value problem (notice that $y_0 > 0$).
We have thus shown that $L$ is a periodic orbit.
Best Answer
Hints: This will guide you through the process and you can figure out how to do this in Matlab.
To find the critical points, you want to simultaneously solve $x' = 0, y' = 0$. You will get two critical points at $$(x,y) = \left(\dfrac{1}{2}, \dfrac{1}{6}\right), (1, 0)$$
You can then determine the types of critical points these are by finding the Jacobian, $J(x, y)$, and evaluating the eigenvalues of the $2x2$ Jacobian. In the phase portrait below, you can see we have a stable and an unstable critical point.
The Phase portrait will show these two critical points and should look something like: