$n=x^{2}+xy+y^{2}\leq |x^{2}|+|xy|+|y^{2}| \leq (|x|+|y|)^{2}$, thus $|x|,|y| \leq \sqrt n$, giving finiteness (and the bound $4n+4\sqrt n$ on the number of solutions). To each pair $(x,y)$ we assign the number $x+y\zeta_{6}\in \mathbb C$, where $\zeta_{6}=\frac{1+i\sqrt{3}}2$ is a primitive $6$th root of unity. This assignment is unique (bijective in fact if $x,y \in \mathbb R$) and $|x+y\zeta_{6}|^{2}=x^{2}+xy+y^{2}$. Thus if $|x+y\zeta_{6}y|^{2}= x^{2}+xy+y^{2}=n$ for some integers, we consider the six numbers $a_{k}+b_{k}\zeta_{k}= \zeta_{6}^{k}(x+\zeta_{6}y)$ for $k=0,1,...,5$. By $|\zeta_{6}|=1$ and $\zeta^{k}=1 \iff 6|k$, we get that all of them are different and have the same square of modulus $n$. Simple calculation shows that all $a_{i}, b_{i}$ are integers (we have just to prove $x,y \in \mathbb{Z}\implies a_{1}, b_{1}\in \mathbb{Z}$, but $a_{1}=-y, b_{1}=x+y$). So we can group the solutions together in couples of six, giving the result. But there is a much stronger result based on factorisation in Eisensteinian integers: The number integer $n>0$ can be written as $x^{2}+xy+y^{2}$ ($x,y \in \mathbb Z$) in exactly \[6\cdot\left(\sum_{\substack{d|n \\ d \equiv 1 \mod 3}}1-\sum_{\substack{d|n \\ d \equiv 2 \mod 3}}1 \right)\] ways. If there is one or more prime divisor $p \equiv 2 \mod 3$ dividing $n$ an odd number of times, then this is $0$, otherwise it equals \[6 \prod_{\substack{p \in \mathbb{P}\\ p \equiv 1 \mod 3}}\left(v_{p}(n)+1\right).\]

The discriminant is $4n-3x^2$ which shows there is a finite number of solutions. if (x, y) is a solution then so is (y-x, -x) and (-y, x-y). (-x, -y) is also a solution and this leads to 6 distinct solutions that can be found for every pair (x, y).