Note that 2023=172⋅7 and let (x,y) be a solution. If we assume that y is not divisible by 17, there exists a number a such that ay≡x(mod17). We then have
0=x2−xy+y2≡y2(a2−a+1)(mod17),
hence a2−a+1≡0(mod17). This implies that a6≡1(mod17), because
a6−1=(a3−1)(a3+1)=(a3−1)(a+1)(a2−a+1).
On the other hand, by Fermat's little theorem, a16≡1(mod17). Then a2≡a2⋅a16≡a18≡(a6)3≡1(mod17), hence a≡±1(mod17). But for these two values of a we get a2−a+1≡2−a≡2∓1≡0(mod17), a contradiction. Hence, y must be divisible by 17.
It follows that x must also be divisible by 17, and we look for solutions (x,y)=(17u,17v). Let N(x,y)=x2−xy+y2 and note that
N(x,y)=4(2x−y)2+3y2.
Then 2023=N(17u,17v)=172N(u,v), iff 7=N(u,v), or equivalently
28=(2u−v)2+3v2.
Because v2≤28/3 and 28 is not a square, to find solutions to this equation we only need to consider v∈{±1,±2,±3}. To a solution with negative v there corresponds one with positive v, obtained by changing the signs of both, u and v. Because 28−3=52, 28−12=42, and 28−27=12, we have these three cases for solutions with positive v
v=1,2u−v=±5,v=2,2u−v=±4,v=3,2u−v=±1.
The possibilities for solutions (u,v) with positive v are therefore
(3,1)(−2,1)(3,2)(−1,2)(2,3)(1,3).
Thus we obtain 12 solutions to the original equation, namely
(x,y)=(17,51),(34,51),(−17,34)
and all the variants (−x,−y), (y,x), (−y,−x) of these.