From your stability analysis for the equilibrium point $(0, 0)$ you can see that there are two eigenvectors that correspond to the two eigienvalues of the linearized system: $$\lambda = 1 \,\, \mapsto \,\, \vec{e}_1 = \begin{bmatrix}1\\0\end{bmatrix} \, \, \text{ and } \,\, \lambda = 0 \,\,\mapsto\,\, \vec{v}_1 = \begin{bmatrix}-1\\1\end{bmatrix}$$
One thing that can be checked is whether the lines defined by these vectors are invariant lines of the system of ODEs. Indeed, the line $y = 0$ and the lines $x+y = 0$ are both invariant for the system. For example, for one to check that $x+y=0$ is an invariant line, one has to check that for any starting point on $x+y=0$ the solution that starts from the given starting point always satisfies the equation $x+y=0$, i.e. it always stays on that line. Same for the horizontal line $y = 0$. Such tests can be performed, for example, by checking whether the normal vectors of each of these lines are orthogonal to the vectors of the vector field, defined by the system of ODEs, for any point on the given lines.
Second Edit.
Definition. Let
\begin{align}
&\frac{dx}{dt} \,=\, f(x,y)\\
&\\
&\frac{dy}{dt} \,=\, g(x,y)
\end{align} be a system of ODEs, where $f$ and $g$ are smooth functions (for simplicity say on all of $\mathbb{R}^2$). Let $\Gamma$ be a smooth curve in $\mathbb{R}^2$, defined as $$\Gamma \,=\, \big\{\,(x,y) \in \mathbb{R}^2 \, :\,\, q(x, y) = 0\,\big\}$$ where $q$ is also a smooth function, and $\frac{\partial q}{\partial x}(x,y)\neq 0$ or $\frac{\partial q}{\partial y}(x,y)\neq 0$ for each point $(x, y) \in \Gamma$. Then, the curve $\Gamma\, :\, \, q(x, y) = 0$ is called invariant curve for the system above whenever, if a solution $\big(x(t), \, y(t)\big)$ of the system of ODEs satisfies the initial condition $q\big(x(t_0), y(t_0)\big) = 0$ for some $t_0 \in \mathbb{R}$, then $q\big(x(t), y(t)\big) = 0$ for all $t \in \mathbb{R}$ for which the solution exists. In other words, an invariant curve is a smooth curve $\Gamma$ such that if a solution trajectory of the system has at least one point in common with the curve $\Gamma$, then the whole solution trajectory lies on the curve $\Gamma$.
Proposition. Let
\begin{align}
&\frac{dx}{dt} \,=\, f(x,y)\\
&\\
&\frac{dy}{dt} \,=\, g(x,y)
\end{align} be a system of ODEs, where $f$ and $g$ are smooth functions (for simplicity say on all of $\mathbb{R}^2$). Let $\Gamma$ be a smooth curve in $\mathbb{R}^2$, defined by the equation $q(x, y) = 0$ where $q$ is also a smooth function and $\frac{\partial q}{\partial x}(x,y)\neq 0$ or $\frac{\partial q}{\partial y}(x,y)\neq 0$ for every point $(x, y) \in \mathbb{R}^2$ such that $q(x, y) = 0$. Then, if for all $(x, y) \in \mathbb{R}^2$ (or in an open subset of $\mathbb{R}^2$ containing the curve $\Gamma$)
$$f(x, y)\,\frac{\partial q}{\partial x}(x,y) \,+\, g(x, y)\,\frac{\partial q}{\partial y}(x,y) \, =\, \lambda(x,y)\, q(x, y)$$ for some continuous function $\lambda(x, y)$, then the curve $\Gamma\, :\, \, q(x, y) = 0$ is invariant curve for the system of ODEs above.
Intuitive idea of the proof: An invariant curve means that if as solution trajectory starts from a point $(x, y)$ on the curve $\Gamma$, then the whole solution trajectory lies on the curve $\Gamma$. Which also means that, for the curve $\Gamma$ to be invariant, the velocity tangent vector to the said solution trajectory on the curve must be tangent to the curve itself. However, a vector at a point $(x, y)$ is tangent to the curve $\Gamma$ if and only if it is orthogonal to the gradient vector $\Big( \,\frac{\partial q}{\partial x}(x,y), \,\frac{\partial q}{\partial y}(x,y)\,\Big)$ of the curve at the point $(x,y)$. Recall that the tangent vectors to the the solution trajectories of the system of ODES above are the vectors defined by the vector field of the system of ODEs
$$f(x,y)\,\frac{\partial}{\partial x} \, +\, g(x,y)\,\frac{\partial}{\partial y}$$ Therefore, the curve $\Gamma$ is invariant whenever the vector field of the system of ODEs is perpendicular to the gradient of the curve $\Gamma$, at the points on the curve. The condition
$$f(x, y)\,\frac{\partial q}{\partial x}(x,y) \,+\, g(x, y)\,\frac{\partial q}{\partial y}(x,y) \, =\, \lambda(x,y)\, q(x, y)$$
reveals exactly that, because for all points on the curve $\Gamma$, $\,\, q(x, y) = 0.$ $\,\,\,\square$
Let us check that this proposition can be applied to the vector field of the system we are studying:
$$\big(\, x + y - x^2\,\big)\,\frac{\partial}{\partial x} \, -\, x\,y\,\frac{\partial}{\partial y}$$
where the curve we are testing for invariance is defined by $q(x, y) \, =\, x + y$. Thus, we need to carry out the following calculation:
\begin{align}
\big(\, x + y - x^2\,\big)\,\frac{\partial}{\partial x}(x + y) \, -\, x\,y\,\frac{\partial}{\partial y}(x + y) \, &=\, \big(\, x + y - x^2\,\big)\cdot 1\, -\, x\,y\cdot 1\\
&=\,\big(\, x + y - x^2\,\big)\, -\, x\,y\\
&=\, x + y - x^2\, -\, x\,y\\
&=\, \, x + y - x\,(x + y)\\
&=\, \, (1 - x)\,(x + y)
\end{align}
which is exactly what we need. Now, by the proposition, the straight line $x + y = 0$ is an invariant line for the system of ODEs.
End of Second Edit.
With this information at hand, and observing that the other equilibrium point is $(1, 0)$, it makes sense to perform a projective transformation of the plane, i.e. a change of variables, that preserves the horizontal line $y = 0$ (it sends it to itself, but transforms its points) and sends the line $x + y = 0$ to the infinite line. A projective transformation that accomplishes such a task is
\begin{align}
& u \, = \, \frac{1}{x+y}\\
& w \, = \, \frac{y}{x+y}
\end{align}
Then, the inverse projective transformation is
\begin{align}
& x \, = \, \frac{1-w}{u}\\
& y \, = \, \frac{w}{u}
\end{align}
When you express the original system of ODEs
\begin{align}
& \frac{dx}{dt} \, = \, \phantom{-}\,x \, +\, y \, -\, x^2\\
&\\
& \frac{dy}{dt} \, = \,-\, x \, y
\end{align}
in these new coordinates $(u, w)$, the result looks like this
\begin{align}
& \frac{du}{dt} \, = \, -\,u \, -\, w \, +\,1\\
&\\
& \frac{dw}{dt} \, = \,-\, w
\end{align}
Indeed,
\begin{align}
& \frac{dx}{dt} \, = \, \frac{d}{dt}\left(\frac{1-w}{u}\right) \,=\, \phantom{-}\,x \, +\, y \, -\, x^2 \, =\, \frac{1-w}{u}
\, + \, \frac{w}{u} \, -\, \left(\frac{1-w}{u}\right)^2\\
&\\
& \frac{dy}{dt} \, = \,\frac{d}{dt}\left(\frac{w}{u}\right) \,=\, -\, x \, y \, =\, \, - \, \frac{1-w}{u} \,\cdot \frac{w}{u} \, =\, \frac{w\,(w-1)}{u^2}
\end{align}
which becomes
\begin{align}
&\frac{d}{dt}\left(\frac{1-w}{u}\right) \, =\, \frac{1}{u}
\, -\, \frac{(1-w)^2}{u^2}\\
&\\
&\frac{d}{dt}\left(\frac{w}{u}\right) \, =\, \frac{w\,(w-1)}{u^2}
\end{align}
When we perform the differentiation on the right
\begin{align}
&\frac{- \,u\,\frac{dw}{dt} - (1-w)\frac{du}{dt}}{u^2} \, =\,
\frac{u - (1-w)^2}{u^2}\\
&\\
&\frac{u\,\frac{dw}{dt} - w\frac{du}{dt}}{u^2} \, =\, \frac{w\,(w-1)}{u^2}
\end{align}
and cancel out the common denominator $u^2$, we get
\begin{align}
&- \,u\,\frac{dw}{dt} - (1-w)\frac{du}{dt} \, =\,
u - (1-w)^2\\
&\\
&u\,\frac{dw}{dt} - w\frac{du}{dt} \, =\, w\,(w-1)
\end{align}
or in matrix form
\begin{align}
\begin{bmatrix}
w-1 & -u\\
-w & u
\end{bmatrix}
\begin{bmatrix}
\frac{du}{dt}\\
\frac{dw}{dt}
\end{bmatrix} \, =\,
\begin{bmatrix}
u - (1-w)^2\\
w\,(w-1)
\end{bmatrix}
\end{align}
so when you inverts the Jacobi matrix on the left, one gets the system
\begin{align}
\begin{bmatrix}
\frac{du}{dt}\\
\frac{dw}{dt}
\end{bmatrix} \, =\,
\begin{bmatrix}
w-1 & -u\\
-w & u
\end{bmatrix}^{-1}\,
\begin{bmatrix}
u - (1-w)^2\\
w\,(w-1)
\end{bmatrix}
\end{align}
To invert the Jacobi matrix, one can simply calculate the determinant first
$$
\det\left(
\begin{bmatrix}
w-1 & -u\\
-w & u
\end{bmatrix}
\right) \, =\, (w-1)\, u \, - \, w\, u \, =\, w\,u - u \, -\, w\,u \,=\, -\,u
$$
and then write
\begin{align}
\begin{bmatrix}
\frac{du}{dt}\\
\frac{dw}{dt}
\end{bmatrix} \, &=\,-\,\frac{1}{u}
\begin{bmatrix}
u & u\\
w & w-1
\end{bmatrix}\,
\begin{bmatrix}
u - (1-w)^2\\
w\,(w-1)
\end{bmatrix}\\
&=\, -\,\frac{1}{u}
\begin{bmatrix}
u\,\Big(\,u - (1-w)^2 + w\,(w-1) \,\Big)\\
w\,\Big(\,u - (1-w)^2\,\Big)\,+\, (w-1) \Big(\,w\,(w-1)\,\Big)
\end{bmatrix}\\
&=\, -\,\frac{1}{u}
\begin{bmatrix}
u\,\Big(\,u - 1 + 2\,w - w^2 + w^2- w \,\Big)\\
u\,w - w\,(1-w)^2\,+\, \,w\,(w-1)^2
\end{bmatrix}\\
&=\, -\,\frac{1}{u}
\begin{bmatrix}
u\,\Big(\,u - 1 + w \,\Big)\\
u\,w
\end{bmatrix}\\
&=\,
\begin{bmatrix}
- u + 1 - w\\
-\, w
\end{bmatrix}
\end{align}
Finally, we end up with the system of ODEs in the new $(u, w)$ coordinates
\begin{align}
& \frac{du}{dt} \, = \, -\,u \, -\, w \, +\,1\\
&\\
& \frac{dw}{dt} \, = \,-\, w
\end{align}
One should have in mind that the line $u = 0$ (the vertical $w$ axis) is the line at infinity of the original system, so it should serve as a separator between the solution trajectories to the left of it and to the right of it. The equilibrium point $(x=1, y=0)$ of the original system is mapped again to $(u=1, w=0)$, so as one can observe, in the new coordinates the system is linear and centered at the equilibrium point $(1, 0)$, where the matrix of the system is
$$
\begin{bmatrix}
-1 & -1\\
\phantom{-}0 & -1
\end{bmatrix}
$$
Edit:
The projective change of variables from $(x, y)$ to $(u, w)$ allows us to observe that when the system is written in the $(u, w)$ coordinates the equilibrium point $(1, 0)$ is asymptotically stable attracting point, with basin of attraction the whole $(u, w)$ plane. However, when one switches back to the original $(x, y)$ coordinates, the vertical line $u = 0$ is transformed to the line at infinity, i.e. the solution trajectories of the original system in $(x, y)$ coordinates do not really get there. Therefore, the basin of attraction of $(1, 0)$ in the $(u, w)$ must be separated by the vertical line $u = 0$. The equilibrium point $(u=1, w=0)$ is in the domain $u > 0$. The domain $u > 0$ is a basin of attraction for $(u=1, w=0)$. Then, the projective transformation from $(u, w)$ back to $(x, y)$ transforms the domain $u > 0$ onto the domain $x + y > 0$: $$u > 0 \,\text{ is the same as }\, u = \frac{1}{x+y} > 0 \, \text{ which is possible if and only if}\, x + y > 0$$
Therefore, the domain $x+y > 0$ is the basin of attraction of the equilibrium point $(x=1, y=0)$. The other half-plane $u < 0$ in the $(u, w)$ coordinates is mapped onto the domain $x + y < 0$. Since $x+y = 0$ is an invariant line for the system, no solution trajectory from $x+y < 0$ can cross it, which means no solutions starting from $x+y < 0$ can ever get to $(1, 0)$ because they cannot cross into $x+y > 0$. Actually, the solutions starting from $x+y < 0$, when viewed in the $(u, w)$ coordinates, travel towards $u=0$, because when extended over $u = 0$ they cross it and go on to wards $(1, 0)$. Hence, in the $(x, y)$ the solutions starting from $x+y < 0$
are attracted to infinity, so they diverge towards infinity and are not attracted to any finite point.
End of Edit.
Finally, you can solve explicitly the linear system in the $(u, w)$ coordinates and then construct the explicit solutions of the original system in $(x, y)$ coordinates, by simply using the coordinate change of variables formulas. Indeed, for the linear system (if I am not mistaken), the solutions are
\begin{align}
&u\,=\, u_0\,e^{-t} - w_0\,t\,e^{-t} + 1\\
&w \,=\, w_0\,e^{-t}
\end{align}
and consequently, the solutions of the original system are
\begin{align}
&x\,=\, \frac{1 - w_0\,e^{-t}}{u_0\,e^{-t} - w_0\,t\,e^{-t} + 1}\\
&y \,=\, \frac{w_0\,e^{-t}}{u_0\,e^{-t} - w_0\,t\,e^{-t} + 1}
\end{align}
Best Answer
In addition to what has been said by @Evgeny and @MrYouMath: the set $$ M=\left\{ (x,y)\in\mathbb R^2 :\; x^2+y^2<2 \right\} $$ is a positively invariant set of the considered system since $\forall (x,y)\in M$ $$ \dot V=-x^4-y^4+x^2y^2(x^2+y^2)<-x^4-y^4+2x^2y^2=-(x^2-y^2)^2\leq 0; $$ it is also a subset (guaranteed estimation) of the domain of attraction.