Here's what I think:
Based on the (presumably correct) given form for the solutions:
$x=e^{{c_1}e^{-t}}+c_2e^{2t}-1, \tag{1}$
$y=c_1e^{-t}, \tag{2}$
we can obtain the stable and unstable manifolds of the system by arguing as follows:
the stable manifold is, by definition, the set of points $p$ such that $\phi_t(p) \to (0, 0)$ as $t \to \infty$, where $\phi_t$ is the flow of the given system. Since for nonzero $c_2$, $c_2e^{2t} \to \text{sign}(c_2) \infty$ as $t \to \infty$, $c_2 = 0$ on the stable manifold. It is thus given by
$x=e^{{c_1}e^{-t}}-1, \tag{3}$
$y=c_1e^{-t}, \tag{4}$
or if you like
$x = e^y - 1, \tag{5}$
or
$y = \ln (x + 1), x > -1. \tag{6}$
Note that from (3), we have $x > -1$ for any combination of finite $c_1$ and $t$, so (6) is well-defined. Likewise, the unstable manifold consists of those $p$ such that
$\phi_t(p) \to (0, 0)$ as $t \to -\infty$. This condition forces $c_1 = 0$ so that the unstable manifold is given by
$x = c_2e^{2t}, \tag{7}$
$y = 0, \tag{8}$
i. e., it is the entire $x$-axis. For a further analysis of this system, see my answer to the OP's other question relating to it, here.
Nota Bene: Not to put too fine a point on it, but it appears to me that there is an erroneous assertion in Evgeny's answer to this question. Though it is true that $y$ is decoupled from $x$, since $\dot y = -y$, the $y$-axis is not invariant under the flow. To see this, consider $\dot x$ at some point on the $y$-axis, $(0, y_0)$. Then $\dot x = 2 - (2 + y_0)e^{y_0} = 2(1 - e^{y_0}) - y_0e^{y_0} <0$ for $y_0 > 0$. This shows that at least some points on the $y$-axis are moved off of it by the flow, so it cannot be invariant. End of Note.
Hope this helps. Cheers,
and as always,
Fiat Lux!!!
First note that the third equation has the solution
$$
z(t)^{-2}-z_0^{-2}=2t\implies z(t)=\frac{z_0}{\sqrt{1+2z_0^2t}}
$$
where now indeed $z(0)=z_0$.
Next, you can be a little creative in solving the second equation by inserting the third to reduce the 4th degree term,
$$
\dot y+2y=2z^2-2z^4=2z^2+2z\dot z=\frac{d}{dt}(z^2)+2z^2.
$$
This means that "luckily" the terms turn out to be that the same differential operator $(D+2)$ is applied to both sides. The solution then is
$$
y(t)-z(t)^2=(y_0-z_0^2)e^{-2t},\\
y(t)=\frac{z_0^2}{1+2z_0^2t}+(y_0-z_0^2)e^{-2t}.
$$
I do not see how there can be any unstable manifold. From points outside the set with $x_0=0$ and $y_0=z_0^2$ the solutions will converge exponentially fast towards points within that set, and then slowly along that set towards the origin. In my mind this was the characteristic of a center manifold?
Best Answer
This is a standard problem in linear dynamical systems. And the easiest way to find stable sets is to find the eigenvalues. According to that, one can categorize whether the origin is a source, sink, saddle etc. The classic reference on this is Hirsch and Smale's Differential equations, dynamical systems and linear algebra. The book has a second edition which includes some chaos theory but I prefer the original first edition.