Finding a Hopf Bifucation with eigenvalues I am trying to show that the following 2D system has a Hopf bifurcation at $\lambda=0$:
\begin{align}                                                                                                                                                                
x' =& y + \lambda x \\                                                                                                                                                       
y' =& -x + \lambda y - x^2y                                                                                                                                                  
\end{align}
I know that I could easily plot the system with a CAS but I wish to analytical methods. So, I took the Jacobian:
\begin{equation}                                                                                                                                                             
  J = \begin{pmatrix} \lambda&1\\-1-2xy&\lambda-x^2\end{pmatrix}                                                                                                             
\end{equation}
My book says I should look at the eigenvalues of the Jacobian and find where the real part of the eigenvalue switches from $-$ to $+$. This would correspond to where the
system changes stability. So I took the $\det(J)$:
\begin{align}                                                                                                                                                                
  \det(J) =& -\lambda x^2 + 2xy + \lambda^2 + 1 = 0                                                                                                                          
\end{align}
I am stuck here with algebra and am not quite sure how to find out where the eigenvalues switch from negative real part to positive real part. I would like to use the
quadratic formula but the $2xy$ term throws me off.
How do I proceed? Thanks for all the help!

Step 1: As jkn stated, the first step is to find the equilibrium such that $\dot x=\dot y=0$. It is easy to obtain the equilibrium  is $(0, 0)$.
Step 2: Compute the eigenvalue around the equilibrium $(0, 0)$.
$$J= \left(\begin{array}{cc}
                                                         \lambda &1\\
                                                       -1&  \lambda\\
                                                       \end{array}\right)$$
Thus the characteristic equation is $$\det (\tilde\lambda I-J)=0$$
The aobve equation admits two eigenvalue $\lambda\pm i$. It is called the imaginary $\omega=1$.
It is obvious that $\lambda:=0$ such that the dynamical system has a pair of pure imaginary roots $\pm i$, moverover 
$$\frac{d\tilde\lambda}{d\lambda}|_{\lambda=0}=\frac{d}{d\lambda}|_{\lambda=0}(\lambda\pm i)=1>0$$
There the dynamical system exhibits a Hopf bifurcation at $\lambda=0$.
Step 3, Compute the Hopf bifurcation. Consider
$$(J-\tilde\lambda I)\left(\begin{array}{cc}
                                                         q_1\\
                                                       q_2\\
                                                       \end{array}\right)=0$$
we get the eigenvector $$\left(\begin{array}{cc}
                                                         q_1\\
                                                       q_2\\
                                                       \end{array}\right)=\left(\begin{array}{cc}
                                                         i\\
                                                       1\\
                                                       \end{array}\right)$$.
On the other hand, consider$$(J^T-\overline{\tilde\lambda} I)\left(\begin{array}{cc}
                                                         p_1\\
                                                       p_2\\
                                                       \end{array}\right)=0$$
we get the eigenvector $$\left(\begin{array}{cc}
                                                         p_1\\
                                                       p_2\\
                                                       \end{array}\right)=\left(\begin{array}{cc}
                                                         i\\
                                                       1\\
                                                       \end{array}\right)$$
It is always possible to normalize $\mathbf p$ with respect to $\mathbf p$:
$$<\mathbf p, \mathbf q>=1, \mbox{ with } <\mathbf p, \mathbf q>=\bar p_1q_1+\bar p_2 q_2 $$
Thus $$\left(\begin{array}{cc}
                                                         p_1\\
                                                       p_2\\
                                                       \end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}
                                                         i\\
                                                       1\\
                                                       \end{array}\right)$$
We denote the nonlinear term $$\mathbf F:=\left(\begin{array}{cc}
                                                         F_1\\
                                                       F_2\\
                                                       \end{array}\right)
=\left(\begin{array}{cc}
                                                         0\\
                                                       -x^2y\\
                                                       \end{array}\right)$$
Introduce the complex variable $$z:=<\left(\begin{array}{cc}
                                                         p_1\\
                                                       p_2\\
                                                       \end{array}\right),
\left(\begin{array}{cc}
                                                         x\\
                                                       y\\
                                                       \end{array}\right)>$$
Then the dynamical system becomes to
$$\dot z=(\lambda+i)z+g(z, \bar z,\lambda) \mbox{ with }g(z, \bar z,\lambda) =<\left(\begin{array}{cc}
                                                         p_1\\
                                                       p_2\\
                                                       \end{array}\right),
\left(\begin{array}{cc}
                                                         F_1\\
                                                       F_2\\
                                                       \end{array}\right)>$$
Some direct computations show
$$g(z, \bar z,\lambda)=\frac{1}{2}(z^3-z^2\bar z-z\bar z^2+\bar z^3):=\frac{g_{30}}{6}z^3+\frac{g_{21}}{2}z^2\bar z+\frac{g_{12}}{2}z\bar z^2+\frac{g_{03}}{6}\bar z^3$$
We recall that the first Lyapunov coefficient is
$$l_1=\frac{1}{2\omega^2}\mbox{Re} (ig_{20}g_{11}+\omega g_{21})$$
Since $\omega=1$, $g_{20}=g_{11}=0$, $g_{21}=-1$, we obtain
$$l_1=-\frac{1}{2}<0$$
Therefore $\lambda=0$ is a super crirical Hopf bifurcation.
