<s>[INST] Simultaneous real solution of $x^3+y^3+1+6xy=0$ & $xy^2+y+x^2=0$ I am trying to solve the following system of non-linear equations in real numbers: 
$x^3+y^3+1+6xy=0$ & $xy^2+y+x^2=0$, with $x,y$ real.
I can only see that $xy\ne 0$. 
I have no clue whether a solution exists or not and how to find any solution. I cannot seem to be able to separate $x,y$ or write it in a good parametric form.
Please help.
 [/INST] The given set of equations has three real solutions:
$$(x,y) =
\left(\frac{1}{y_2},y_1\right),
\left(\frac{1}{y_3},y_2\right),
\left(\frac{1}{y_1},y_3\right)
$$
where $y_1 = 2\cos\frac{\pi}{9}$, $y_2 = 2\cos\frac{7\pi}{9}$ and $y_3 = 2\cos\frac{13\pi}{9}$. 
To derive them, we will exploit a symmetry hidden underneath the set of equations. 
Introduce variables $(X,Y,Z)$ and let $(x,y) = \left(\frac{X}{Z},\frac{Y}{Z}\right)$, 
we have
$$\begin{align}
x^3 + y^3 + 1 + 6xy = 0 
& \iff X^3 + Y^3 + Z^3 + 6XYZ = 0 \tag{*1a}\\
xy^2 + y + z^2 = 0 
&\iff XY^2 + YZ^2 + X^2Z = 0\tag{*1b}
\end{align}
$$
Let $u = \frac{Y}{Z}, v = \frac{Z}{X}, w = \frac{X}{Y}$, we have $uvw = 1$. 
Divide RHS(*1a) by $XYZ$, we obtain
$$\frac{w}{v} + \frac{u}{w} + \frac{v}{u} +  6 = \frac{X^2}{YZ} + \frac{Y^2}{XZ} + \frac{Z^2}{XY} + 6 = 0$$
Multiply both sides by $uvw = 1$, we obtain
$$w^2 u + u^2v + v^2w + 6 = 0 \tag{*2}$$
Similarly, divide RHS(*1b) by $XYZ$ give us
$$u + v + w = \frac{Y}{Z} + \frac{Z}{X} + \frac{X}{Y} = 0$$
Let $\omega = e^{\frac{2\pi}{3}i}$ be the cubic root of unity. Since $u+v+w = 0$, we can find
two complex numbers $a, b$ such that
$$\begin{cases}
u &= a + b\\
v &= a\omega + b\omega^2\\
w &= a\omega^2 + b\omega
\end{cases}\tag{*3}$$
In terms of $a,b$, the condition $uvw = 1$ reduces to
$$a^3 + b^3 = (a+b)(a\omega + b\omega^2)(a\omega^2+b\omega) = uvw = 1$$
Substitute the parameterization (*3) into (*2) and simplify, we obtain
$$3a^3\omega + 3b^3\omega^2 + 6 = 0 \iff a^3\omega + b^3\omega^2 = - 2$$
Solving this last two equations give us
$$(a^3, b^3 )  = ( -\omega^2, -\omega ) = \left( e^{\frac{\pi}{3}i}, e^{-\frac{\pi}{3}i}\right)$$
We are looking for real solutions for $(x,y)$, we need $u = a + b$ and $v = a\omega + b\bar{\omega}$ to be real. This forces $b = \bar{a}$ and there are three possibilities:
$$a = \bar{b} = e^{\frac{\pi}{9}i}, e^{\frac{7\pi}{9}i}, e^{\frac{13\pi}{9}i}
\implies (u,v,w) = (y_1,y_2,y_3), (y_2,y_3,y_1), (y_3,y_1,y_2)$$
Notice $(x,y) = (v^{-1},u)$, the three real solutions of the problem follow.
</s>