Hersh's answer solves the problem. I want to point out, though, a general recipe you can follow for solving these kinds of problems:
- Use trigonometric identities (in particular angle-sum and multiple-angle identities) to transform your equations into polynomials in sines and cosines of the unknown variables.
In your case you will get:
\begin{align*}
\cos x + \cos x \cos y - \sin x \sin y &= 0\\
\cos y + \cos x \cos y - \sin x \sin y &= 0.
\end{align*}
In general, one case when this transformation is always possible is when your trig function arguments are linear polynomials in the unknown variables with integer coefficients.
- Use the substitution $c_x = \cos x, s_x = \sin x$, etc. to replace all trig functions with unknown real variables. Add an equation of the form $c_x^2+s_x^2=1$ for each variable. This new equation guarantees that if you can solve for $c_x$ and $s_x$, you will be able to find an $x$ whose cosine and sine are equal to the variable values.
In your case you will get
\begin{align*}
c_x + c_x c_y - s_x s_y &= 0\\
c_y + c_x c_y - s_x s_y &= 0\\
c_x^2+s_x^2 &= 1\\
c_y^2+s_y^2 &= 1.
\end{align*}
- Use standard algebra techniques (elimination, substitution, etc.) to solve the system of polynomial equations. This is the one part of the process that's not pure recipe and may require some insights to avoid nasty calculation.
In your case, you can follow essentially Hersh's approach and subtract the first two equations to get
$$c_x = c_y.$$
Now the third and fourth equations give you
$$s_x = \pm\sqrt{1-c_x^2},\qquad s_y = \pm \sqrt{1-c_x^2}.$$
Now plug these into the first equation. When $s_x$ and $s_y$ have the same sign, you get
$$c_x + c_x^2 - (1-c_x^2) = 0$$
a quadratic with roots $c_x \in \left\{-1, \frac{1}{2}\right\}.$ These give you three solutions:
$$(c_x, s_x, c_y, s_y) = \left(-1, 0, -1, 0\right), \left(\frac{1}{2}, -\frac{\sqrt{3}}{2}, \frac{1}{2}, -\frac{\sqrt{3}}{2}\right), \left(\frac{1}{2}, \frac{\sqrt{3}}{2}, \frac{1}{2}, \frac{\sqrt{3}}{2}\right).$$
When $s_x$ and $s_y$ have opposite signs, you get $c_x + c_x^2 + (1-c_x^2) = 0$, which gives you $(-1, 0, -1, 0)$ again and no new solutions.
- Recover the solution for the original variables by computing $x = \arctan(s_x, c_x) + 2\pi n_x$. (Here we need to use the "two-argument arctangent" that correctly accounts for the quadrants of the inputs.)
For your problem we get (up to multiples of $2\pi$):
$$(x,y) = (\pi,\pi), \left(-\frac{\pi}{3}, -\frac{\pi}{3}\right), \left(\frac{\pi}{3}, \frac{\pi}{3}\right).$$
Finally for your particular problem you want the solution in $[0,2\pi)$ so of course you'll want $\frac{5\pi}{3}$ instead of $-\frac{\pi}{3}.$