Diophantine equation $x^2 + xy + y^2 = \left({{x+y}\over{3}} + 1\right)^3$. Solve in integers the equation$$x^2 + xy + y^2 = \left({{x+y}\over3} + 1\right)^3.$$

Setting $x+y=3t$ and $xy=s$, we obtain that
$$9t^2-s = (t+1)^3 \implies s = -t^3+6t^2-3t-1$$
$x$ and $y$ satisfying the quadratic $a^2 -3ta + s =0$. This means $9t^2-4s = k^2$. Eliminating $s$, we obtain
$$4t^3-15t^2+12t+4 = k^2 \implies 64t^3 - 240t^2 + 192t + 64 = (4k)^2$$
$$(4t-5)^3 - 108t + 189 = (8k)^2 \implies (4t-5)^3 - 27(4t-5) + 54 = (4k)^2$$
Hence, it boils down to solving for integer points on the elliptic curve
$$Y^2 = X^3 - 27X+54 = (X-3)^2(X+6)$$
Hence, we need $X+6=m^2$, i.e., $4t+1=m^2 \implies t = \dfrac{m^2-1}4$. This means $m$ has to be odd. Choosing $m=2p+1$, we obtain
$$t = \dfrac{4p^2+4p}4 = p^2 + p$$
This gives us that
$$Y = m(m^2-9) = (2p+1)(4p^2+4p-8) \implies k=(2p+1)(p+2)(p-1)$$
This gives us $x,y = \dfrac{3t \pm k}2 = t + \dfrac{t \pm k}2$.
Hence, to summarize,


*

*Pick any $p \in \mathbb{Z}$.

*Set $t=p(p+1)$ and $k = (2p+1)(p+2)(p-1)$

*This gives us $\boxed{\color{blue}{x,y = t + \dfrac{t \pm k}2}}$


If we only want unordered pairs, it suffice to vary $p$ just over $\mathbb{N}$. The first few values by varying $p$ from $0$ to $5$ are
\begin{array}{|c|c|c|}
\hline
p & x & y\\
\hline
0 & -1 & 1\\
1 & 3 & 3\\
2 & 19 & -1\\
3 & 53 & -17\\
4 & 111 & -51\\
5 & 199 & -109\\
\vdots & \vdots & \vdots\\
\hline
\end{array}
