Hint.
A good stream plot can help a lot!
In red the equilibrium points for $\epsilon = 0.1$.
The corresponding MATHEMATICA script.
epsilon = 0.1;
str = StreamPlot[{x (x - 0.5) (x + 0.5) + epsilon y,
y (y - 0.5) (y + 0.5) + epsilon x}, {x, -1, 1}, {y, -1, 1},
MeshFunctions -> Function[{x, y, vx, vy, n}, n], Mesh -> 55];
sols = Quiet[Solve[{x (x - 0.5) (x + 0.5) + epsilon y == 0,
y (y - 0.5) (y + 0.5) + epsilon x == 0}, {x, y}]];
equil = Table[Graphics[{Red, Disk[({x, y} /. sols[[k]]), 0.02]}], {k, 1,
Length[sols]}];
Show[str, equil]
Your original system can be rewritten as:
$$\begin{bmatrix} \dot{y}_1 \\ \dot{y}_2 \end{bmatrix}=\begin{bmatrix} -2 & 1 \\ -1 & 0 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}$$
Where $y_1=y_1(t)$, $y_2=y_2(t)$ are the functions to be found.
Let's consider a general system of linear ODEs with constant coefficients:
$$\dot{\vec{y}}=\hat{A} \vec{y}$$
Indeed, a partial solution to such a system can be expressed through the eigenvalues and eigenvectors of a square matrix $\hat{A}$ because:
$$\hat{A} \vec{a}_j=v_j \vec{a}_j$$
$$\vec{y}_j=\vec{a}_j e^{v_j t}$$
Where $\vec{a}_j$ and $v_j$ are eigenvectors and eigenvalues respectively. $j=1,...,n$ where $n$
Using this, we can write:
$$\dot{\vec{y}}_j=v_j\vec{a}_j e^{v_j t}=\hat{A} \vec{a}_j e^{v_j t}=\hat{A}\vec{y}_j$$
We can write the general solution as a combination of partial solutions for all eigenvalues and eigenvectors:
$$\vec{y}(t)=\sum_j^n C_j \vec{a}_j e^{v_j t}$$
For some arbitrary constants $C_j$.
Now eigenvalues for small $n$ are indeed found by solving the equation:
$$\det (\hat{A}-v\hat{I})=0$$
In this case the equation becomes:
$$(-2-v)(-v)+1=0 \\ v^2+2v+1=0 \\ (v+1)^2=0 \\ v=-1$$
We get a degenerate case, when the two eigenvalues are the same (repeated eigenvalues).
Because of that, for the general solution we need to take the following two vectors:
$$\vec{y}_1=\vec{a} e^{vt} \\ \vec{y}_2=(\vec{b}+\vec{a} t)e^{vt}$$
Where $\vec{b}$ is the solution of:
$$(\vec{A}-v \vec{I}) \vec{b}=\vec{a}$$
Why can we do that I leave for you to find in your lectures/textbooks.
This leaves us to find the eigenvector $\vec{a}=(a_1,a_2)^T$ for $v=-1$. We need to solve:
$$-2a_1+a_2=-a_1 \\ -a_1=-a_2$$
We get $a_1=c$, $a_2=c$ for some arbitrary constant $c$. Since we are searching for a partial solution we can take $c=1$.
The vector $\vec{b}$ is found from:
$$(-2+1) b_1+b_2=1 \\ -b_1+b_2=1$$
Again, we get multiple solutions, so I think we can choose whatever we want, like:
$$b_1=0 \\ b_2=1$$
The general solution can be expressed as:
$$\vec{y}(t)=C_1 \begin{bmatrix} 1 \\ 1 \end{bmatrix} e^{-t}+C_2 \left(\begin{bmatrix} 0 \\ 1 \end{bmatrix}+\begin{bmatrix} 1 \\ 1 \end{bmatrix} t \right) e^{-t}$$
For arbitrary constants $C_1,C_2$. We can check if it works by substitution:
$$\dot{\vec{y}}(t)=(C_2-C_1) \begin{bmatrix} 1 \\ 1 \end{bmatrix} e^{-t}-C_2 \left(\begin{bmatrix} 0 \\ 1 \end{bmatrix}+\begin{bmatrix} 1 \\ 1 \end{bmatrix} t \right) e^{-t}$$
$$\hat{A} \vec{y}=-C_1 \begin{bmatrix} 1 \\ 1 \end{bmatrix} e^{-t}+C_2 \left(\begin{bmatrix} 1 \\ 0 \end{bmatrix}-\begin{bmatrix} 1 \\ 1 \end{bmatrix} t \right) e^{-t}$$
The two expressions after simplification are equal, so our solution is correct.
Best Answer
If you write your solution now as $$\vec{x}(t)=\left(C_1\begin{bmatrix} 1 \\ 1 \end{bmatrix}+C_2\begin{bmatrix} 0 \\ 1\end{bmatrix}\right)e^{-t}+C_2\begin{bmatrix}1 \\ 1\end{bmatrix}te^{-t}$$
it may be a little easier to see what's going on. As $t\to\pm\infty$, the equation is dominated by the second term, so $$\vec{x}(t)\approx C_2\begin{bmatrix} 1 \\ 1\end{bmatrix}te^{-t}\;\;\;\;\;\;\;\;\;\;\;t\to\pm\infty$$
As $t\to\infty$, the whole thing decays to $0$, so it makes sense that all solution curves begin to look like they fall along this eigenvector as they decay to $0$. Applying the same logic as $t\to-\infty$, it would seem that solutions grow along this vector as well, but this is only true on a macroscopic scale, as the solutions are growing exponentially in both directions. The first term, no longer decaying to $0$, shifts the solution away from the eigenvector and is growing large, but not as large as the second term. This means that when zoomed in, solutions look like they are leaving the eigenvector, but zoomed out solutions look like they are along the eigenvector (that is, until you let time grow large enough, in which you will see it move away slowly).
The last thing to consider is the direction of the curve. Suppose that $\vec{x}$ is in the first quadrant as $t\to\infty$. This means that $C_2$ is positive, so if $t\to-\infty$, then the sign on the second term is now negative, and the solution turns around to enter the third quadrant. This means that solutions sort of "spin" by turning $180^\circ$ around while still growing away from the origin.
The vector $\begin{bmatrix} 0 \\ 1\end{bmatrix}$ tells you which way the solution spins. If $C_2$ is positive, then that means that at $t=0$, the solution is above (in the xy-plane) of the line through the vector $\begin{bmatrix} 1 \\ 1\end{bmatrix}$ and the origin, so the solution must remain on this side of said line. The dominating term is in the first quadrant as $t\to\infty$, so the solution must spin clockwise to decay this way. It turns out that if solutions spin clockwise on one side of the eigenvector, they do the same on the other side as well; the same is true if they spin counterclockwise.
This type of equilibrium is called a Degenerate Node, in case you were curious. Here is an example of what one might look like.