Region of attraction of simple ODE with perturbation

Question

There are a few nice discussions about ROA covering a few subtopics:

1. Region of attraction of : $x'=-y-x^3,y'=x-y^3$ via Lyapunov Function
2. Region of attraction and stability via liapunov's function (No2)
3. Basin of attractions with continuity and discontinuity
4. Strength of attraction of fixed points
5. Can the basin of attraction be multiple intervals?

Now, consider the following simplest example:

$$\dot{x} = x(x-1)(x+1)$$ $[-1,1]$ is the ROA.

Now consider the two dimensional case:

\begin{equation} \begin{aligned} &\dot{x} = x(x-1)(x+1)\\ &\dot{y} = y(y-1)(y+1) \end{aligned} \end{equation}Obviously, ROA is a square. However, if I consider the following coupled ODE:

\begin{equation} \begin{aligned} &\dot{x} = x(x-1)(x+1) + \epsilon (y-x)\\ &\dot{y} = y(y-1)(y+1) + \epsilon (x-y) \end{aligned} \end{equation} where $\epsilon$ is a very small number. Or

\begin{equation} \begin{aligned} &\dot{x} = x(x-1)(x+1) + \epsilon (-y+x)\\ &\dot{y} = y(y-1)(y+1) + \epsilon (-x+y) \end{aligned} \end{equation} Then I have the following ROAs: (blue line-case three, black line-case two, red line-case one)

My questions are:

1. There are two different tilt directions for case two and three. I know this is because of the slope of the coupling term (for case two, the slope of $x$ and $y$ in the coupling terms are $-1$). But how could I analyze this formally?

2. Is it a good way to analyze 1. by perturbation method (I think this should be a regular perturbation method)? (observe the sign of the leading order term of the solution obtained from perturbation method?) and how could I proceed it for the coupling term?

Note: It is simple to check that if you just use the linearization method to find the Jacobian matrix (w.r.t the point $(0,0)$), the ROA will be the whole $\mathbb{R}^2$, which is not correct.