Since the method below is probably correct, and correctness is potentially irrelevant to my ability to do what I want to learn. Assume below is correct. $\def\b{\begin{pmatrix}}\def\e{\end{pmatrix}}\def\1{y_1}\def\2{y_2}\def\R{\Rightarrow}\def\l{\lambda}\def\le{\left}\def\ri{\right}$
$\dot{y}= \b 3\1+\1\2 \\ \1+\2 -\2^2 \e$
This occurs when $\dot{y}=0$.
Converting $\dot{y}$ to a nicer form:
$\dot{y}= Ay+u$, where $A = \b 3 & \1 \\ 1&1\e, u = \b 0\\-\2^2 \e$
No such nicer form apparently.
So now I just 'go it': $$3\1+\1\2=0\;\;(1)$$ $$\1+\2-\2^2=0\;\;(2)$$
$$(1)\R \1(3+\2)=0, \1=0 \text{ or } \2 = -3$$ $$(2)\R \1=\2(\2-1)$$
Hence we have: $(0,0),(0,1),(12,-3)$
Now I believe I take the eigenvectors of $A$ given these $\1$.
$\1=0,\;\;\det\le| A \ri| = \det\le|\b 3-\l &0\\1&1-\l\e\ri|=(3-\l)(1-\l)\R\l=1,3$
Which is an unstable improper node.
$\1=12,\;\;\det\le| A \ri| = \det\le|\b 3-\l &12\\1&1-\l\e\ri|=(3-\l)(1-\l) -12\R 3 -3\l-\l+\l^2=0\R\l^2-4\l+3\R\l=-3,-1$
This is a stable improper node at both $(0,0)$ and $(0,1)$.
How do I derive the linear system about each critical point: $(0,0),(0,1),(12,-3)$
The linearization of a system $\dot y=F(y)$ around an equilibrium $y_*$ is $$\dot y = F'(y_*)(y-y_*)\tag1$$ where $F'$ is a square matrix (Jacobian). This is based on the fact that $F(y)\approx F'(y_*)(y-y_*)$ neat $y_*$. It is convenient to let $z=y-y_*$, thus rewriting (1) as $$\dot z = F'(y_*)z\tag2$$ which is a linear system.
In your case $$F(y) = \begin{pmatrix} 3y_1+y_1y_2 \\ y_1+y_2 -y_2^2 \end{pmatrix}$$ hence $$F'(y) = \begin{pmatrix} 3+y_2 & y_1 \\ 1 & 1- 2y_2 \end{pmatrix}$$ Plug in the equilibrium point and you get linear system, the stability of which can be studied by looking at the eigenvalues of $F'(y_*)$.
No, you take the eigenvalues of $F'(y_*)$. Only at the zero equilibrium these are the same thing. At other equilibria, $-y_2^2$ has nonzero derivative and therefore contributes to the linearized system.