please solve quadratic simultaneous equations. $x,y,z$ are real variables. $l,m,n$ are real positive constants.
Solve for real $x,y,z$:
$$\begin{align}
x^2+xy+y^2 &= l\\
y^2+yz+z^2 &= m\\
z^2+zx+x^2 &= n
\end{align}
$$
This problem has a relevance in Electrical Engineering.
 A: $$x^2+xy+y^2=l~~(1),~ y^2+yz+z^2=m ~~~(2),~ z^2+zx+x^2=n~~~(3)$$
Subtract (1) from (2) and (1) from (3) to get
$$(x-y)(x+y+z)=m-l~~~(4),~ (x-z)(x+y+z)=n-l~~~(5)$$
Let $(x+y+z)=w$, then (4) and (5) give
$$(x-y)=\frac{m-l}{w}~~~(6), ~(x-z)=\frac{n-l}{w}~~~(7)$$
From these two we get $$[3x-(x+y+z)]w=m+n-2l \implies x=\frac{w^2+m+n-2l}{3w}~~~(8).$$
Similarly, $$y=\frac{w^2+l+b-2m}{3w}, ~~z=\frac{w^2+l+m-2n}{3w}~~~(9)$$
Putting these in Eq. (4), we get
$$w^4-(l+m+n)w^2+(l^2+m^2+n^2)-lm-mn-nl=0$$
This gives 
$$w^2=\frac{(l+m+n)\pm\sqrt{(l+m+n)^2-4[l^2+m^2+n^2-lm-mn-nl]}}{2}$$
$$w^2=\frac{(l+m+n) \pm \sqrt{3[2(lm+mn+nl)-(l^2+m^2+n^2)}}{2}~~~~(10)$$
Lastly, Eqs. (8), (9) and (10) give the values of $x,y,z$ in terms of$l,m,n.$. The roots will be real if $$l^2+m^2+n^2 \le 2(lm+mn+nl)$$
A: Edit : I had written a first solution, now placed at the bottom of this answer.
In fact, there is a more thorough way to work on this issue, with a geometric view that helps to understand differently the kind of connection between the different constraints. Here it is.
First of all, in order to avoid confusions, let us write the system with different letters :
$$\begin{cases}
a^2+ab+b^2 &= l & \ (i)\\
b^2+bc+c^2 &= m & (ii)\\
c^2+ca+a^2 &= n & (iii)
\end{cases}\tag{1}$$
Remark : if $(a,b,c)$ is a solution $(-a,-b,-c)$ is a solution as well.
Consider at first the following figure :

Fig. 1 : The 3 ellipses with equations $x^2+xy+y^2=l,m,n$ in the case $l=7,\color{blue}{m=28},\color{red}{n=21}$ and the hexagon (with  horizontal, vertical or $-45°$ oriented sides) "binding" solutions $(a,b),(b,c),(c,a)$ of the 1st, 2nd, 3rd. equation resp. One must understand that the hexagon must have two of its diagonals horizontal and vertical : in particular, if one moves up a little point $(a,b)$ on the black ellipse, the horizontal diagonal will no longer be acceptable.
_
What are the degrees of freedom we have for determining 3 numbers $a,b,c$ such that 
$$\begin{cases}
(a,b) \in \text{ellipse} (E_l) \text{with equation} \ x^2+xy+y^2 &= l\\
(b,c) \in \text{ellipse} (E_m) \text{with equation} \ x^2+xy+y^2 &= m\\
(c,a) \in \text{ellipse} (E_n) \text{with equation} \  x^2+xy+y^2 &= n
\end{cases} \ ?\tag{2}$$
Let us fix a value of $a$, the associated values of $b$ are solutions to quadratic equation :
$$a^2+ax+x^2=l$$
giving 
$$b=\frac12 (-a \pm \sqrt{4l-3a^2})\tag{3}$$
We can consider, among the two solutions, the positive one corresponding to the first vertex of the hexagon, i.e., point $(a,b) \in (E_l)$. Now,  proceeding on our hexagon in the positive orientation, we reach point  $(c,a)$. How do we find $c$ ? Simply by looking for points of the form $(x,a)$ on curve $(E_n)$ which amounts, in the same way as before, to solve a quadratic equation, here : 
$$x^2+ax+a^2=n$$
giving 
$$c=\frac12 (-a \pm \sqrt{4n-3a^2})\tag{4}$$
Now, comes the crucial moment. We proceed to point $(b,c)$. We have to express that $(b,c) \in (E_m)$, i.e., equation (ii) has to be fulfilled.
Plugging into relationship (ii) the expressions (3) and (4) of $b$ and $c$, we get an equation in the single variable $a$. This equation gives all possible values for $a$, therefore, by (3) and (4), all corresponding values of $b$ and $c$ as wel (in fact, there are $4$ such equations, taking into account the two $\pm$ signs).
This can be done in a numerical way (this is what I have done for testing the method).
One can consider also expand this equation in a smart way by eliminating all square roots by successive squarings in order finally to get a polynomial equation. I haven't succeeded till now.

The previous solution:
We are going to process by necessary conditions.
Multiplying the 1st, 2nd, 3rd equations by $(x-y)$,$(y-z)$,$(z-x)$ resp. , one gets the orthogonality of vector $(l,m,n)$ with vector $V=((x-y),(y-z),(z-x))$. 
Let (P) be the plane orthogonal to vector $(l,m,n)$. A basis of (P) (among others !) is 
$$V_1=(m,-l,0) \ \ \text{and} \ \ V_2=(0,-n,m).$$ 
Let us decompose $V$ onto this basis:
$$V=aV_1+bV_2 \ \ \iff \ \ \begin{cases}x&-&y&&&=&am\\
&&y&-&z&=&-al-bn\\-x&&&+&z&=& \ \ \ \ \ \ \ \ \ \ \ bm\end{cases} \tag{1}$$
adding the 3 equations, we get 
$$0=a(l-m)+b(m-n)$$
which means that $b=a\dfrac{m-l}{m-n}$
Plugging this expression of $b$ into the 3 equations (1), we obtain a parametric solution $(x=x(a),y=y(a),z=z(a))$. 
Now, plug these expressions in the initial equations : it will give 3 quadratic equations say $Q_1=0, \ Q_2=0, \ Q_3=0$ in variable $a$ and it remains to express that these equations have a common real root (see Addendum below)
Up to you for the final computations...
**Addendum : **2 quadratic equations have a common root if their resultant is $0$ ; thus, having a common root will be expressed by 
$$Res(Q_1,Q_2)=0, \ Res(Q_2,Q_3)=0, \ Res(Q_3,Q_1)=0.$$
