Lyapunov first method and Lyapunov function Find if the solutions of following systems are stable, asymptotically stable, or not stable:
a) $\left\{ \begin{array}{ll}
x'=x^3-y \\
y'=x+y^3 \\
\end{array} \right.$
b) $\left\{ \begin{array}{ll}
x'=y-x+xy \\
y'=x-y-x^2-y^3 \\
\end{array} \right.$
c) $\left\{ \begin{array}{ll}
x'=x^5+y^3 \\
y'=x^3-y^5 \\
\end{array} \right.$
d) $\left\{ \begin{array}{ll}
x'=xy-x^3+y^3 \\
y'=x^2-y^3 \\
\end{array} \right.$
I need to use a Lyapunov first method, and if it won't work, then find a Lyapunov function. 
I'm completely new to this (started reading and learning about this today), and can't find full explanation how the first method works. I've seen other examples but not these. Finding Lyapunov function also seems pretty impossible to me. I know how to find eigenvalues and eigenvectors but then what? 
How to build a matrix, when there are multiple x or y terms like in b) or d)?
In a) I've tried obviously $V=x^2+y^2$ and after differentiation I've got $V'=x^4+y^4$ and I think that means, that it's not stable, but the answer says that's stable. 
Will be grateful for any help.
 A: As my explanation is much too long for a comment, I am giving it as an answer.
When you look at a picture of the phase space

you see that it is rather not an example of a stable equilibrium: everywhere the field vectors are pointing outside.
Let's try linearizing.  The Jacobi matrix at $(x, y)$ has the form
$$
J(x, y) = \begin{bmatrix}
3x^2 & -1 
\\
1 & 3y^2
\end{bmatrix}
,
$$
so 
$$
J(0, 0) = \begin{bmatrix}
0 & -1 
\\
1 & 0
\end{bmatrix}
$$
has two eigenvalues, $\pm i$, with zero real parts.  The linearization technique is, regarding stability, inconclusive.
Turn to your Lyapunov function, $V(x, y) = x^2 + y^2$.  When you differentiate $V$ along solutions of the ODE system, you will get
$$
V'(x,y) = 2x \cdot (x^3 - y) + 2y \cdot (x + y^3) = {\color{red}2} (x^4 + y^4),
$$
which is positive everywhere.  What does that mean? Since $V(x,y)$ is the square of the distance from the equilibrium, its derivative along the solutions being positive just means that the distances increase.
A modification of the above argument is as follows.  Consider polar coordinates.  We have 
$$
\begin{cases}
r' \cos{\varphi} - r \sin{\varphi} = r^3 \cos^3{\varphi} - r \sin{\varphi}
\\
r' \sin{\varphi} + r \cos{\varphi} = r \cos{\varphi} + r \sin^3{\varphi}
\end{cases}
$$
which gives
$$
r' = r^3 (\cos^4{\varphi} + \sin^4{\varphi}).
$$
Observe that $\cos^4{\varphi} + \sin^4{\varphi} \ge \frac{1}{2}$, so the function $t \mapsto r(t)$ is increasing.  Indeed, one can show that 
$$
r(t) \ge \frac{r(0)}{\sqrt{1 - (r(0))^2 t}}
$$
as long as a solution exists.  It follows then that whatever small $\delta \in (0, 1)$ we take, if $\sqrt{(x(0)^2+y(0)^2} = \delta$ then there exists $T \in (0, \frac{1}{\delta^2}-1]$ such that $\sqrt{(x(T)^2+y(T)^2} = 1$.  When you compare the quantifiers, you will see that this is even more than instability.
So, if your textbook states that $(0,0)$ is stable, then either that is a misprint, or else ...
