We have the system \begin{align} \begin{cases} \dot{x} = (1-x)^2-y^2\\ \dot{y} = \epsilon^{xy}-1 \end{cases} \end{align} and I must plot the phase graph.
To do this I already know that I must linearize in the first order, since for Taylor degree $\geq2 \ $, $\Delta x \to 0$. However I am stuck and cannot even start.
I have in my mind that I must be able to set the system in such a way such that \begin{align} \begin{cases} y=\text{something}\\ \dot{y}= \text{other something} \end{cases} \end{align} so that then I can run \begin{align} F(x+h)=F(x_0)+h\cdot \dfrac{dF}{dx}|_{x=x_0}. \end{align}
The critical points of the starting system are $(0,1), (0,-1), (1,0)$.
If we have the system
$$\dot{\mathbf{x}}=\mathbf{F}(\mathbf{x})$$
we can linearize this system in point $\mathbf{x}_0$ by using the Jacobian evaluatd at the point $\mathbf{x}_0$ the following method.
$$\dfrac{d}{dt}\left[\mathbf{x}-\mathbf{x}_0\right]=J_F|_{\mathbf{x}=\mathbf{x}_0}\left[\mathbf{x}-\mathbf{x}_0\right]$$
For your system this works as follows for your equilibrium point $\mathbf{x}_0=[x_0,y_0]^T$:
$$\dfrac{d}{dt}\left[x-x_0\right]= \left.\dfrac{\partial}{\partial x}\left[(1-x^2) -y^2 \right]\right|_{x=x_0, y=y_0}\left[x-x_0\right] + \left.\dfrac{\partial}{\partial y}\left[(1-x^2) -y^2 \right]\right|_{x=x_0, y=y_0}\left[y-y_0\right]$$ $$\dfrac{d}{dt}\left[y-y_0\right]= \left.\dfrac{\partial}{\partial x}\left[\epsilon^{xy} -1\right]\right|_{x=x_0, y=y_0}\left[x-x_0\right] + \left.\dfrac{\partial}{\partial y}\left[\epsilon^{xy} -1 \right]\right|_{x=x_0, y=y_0}\left[y-y_0\right]$$
In order to differentiate $\epsilon^{xy}$ you can rewrite it to $\epsilon^{xy}=\exp\left[xy\ln\epsilon \right]$. Executing the partial derivatives and inserting $\mathbf{x}_0=[x_0,y_0]^T$ will lead to
$$\dfrac{d}{dt}\left[x-x_0\right]= -2(1-x_0)\left[x-x_0\right] -2y_0\left[y-y_0\right]$$ $$\dfrac{d}{dt}\left[y-y_0\right]= y_0\epsilon^{x_0y_0}\ln\epsilon^{}\left[x-x_0\right] +x_0\epsilon^{x_0y_0}\ln\epsilon\left[y-y_0\right]$$
We often replace $x-x_0=\Delta x$ and $y-y_0=\Delta y$. Hence, we can rewrite the equations as:
$$\Delta \dot{x}= -2(1-x_0)\Delta x -2y_0\Delta y$$ $$\Delta \dot{y}= y_0\epsilon^{x_0y_0}\ln\epsilon^{}\Delta x +x_0\epsilon^{x_0y_0}\ln\epsilon\Delta y$$
or in matrix form
$$ \dfrac{d}{dt}\Delta\mathbf{x} = \begin{bmatrix} -2(1-x_0) & -2y_0\\ y_0\epsilon^{x_0y_0}\ln\epsilon^{} & x_0\epsilon^{x_0y_0}\ln\epsilon \end{bmatrix}\Delta\mathbf{x}. $$