1. Première méthode :
    1. En développant et en complétant les carrés, on a : \begin{align*} x^2(x^2+7) &= y(2x+y) \\ x^4+7x^2 &= 2xy+y^2 \\ x^4+8x^2+16 &= x^2+2xy+y^2+16 \\ (x^2+4)^2 &= (x+y)^2+16 \qquad (2) \end{align*}

    2. L'équation $(2)$ s'écrit sous forme d'une différence de carrés : \[ ((x^2+4)-(x+y))((x^2+4)+(x+y))=16 \] Les remarques suivantes s'imposent :
      • Les deux facteurs ont la même parité (leur somme est $2(x^2+4)$, qui est paire). Puisque leur produit ($16$) est pair, ils sont nécessairement tous les deux pairs.
      • On a l'inégalité large : $((x^2+4)-(x+y)) \leq ((x^2+4)+(x+y))$ car $x+y \ge 0$.
      Séparons l'étude en deux cas :
      • Cas $x+y=0$ : Puisque $x, y \in \mathbb{N}$, cela implique $x=0$ et $y=0$. On vérifie que $(0^2+4)^2 = 0^2 + 16$, soit $16=16$. Le couple $(0,0)$ est donc une solution.
      • Cas $x+y>0$ : L'inégalité devient stricte : $((x^2+4)-(x+y)) < ((x^2+4)+(x+y))$. Le seul système possible parmi les diviseurs pairs de $16$ est : \begin{align*} \begin{cases} (x^2+4)-(x+y) &= 2 \\ (x^2+4)+(x+y) &= 8 \end{cases} &\implies \begin{cases} 2(x^2+4) &= 10 \\ 2(x+y) &= 6 \end{cases} \\ &\implies \begin{cases} x^2+4 &= 5 \\ x+y &= 3 \end{cases} \\ &\implies \begin{cases} x^2 &= 1 \\ x+y &= 3 \end{cases} \end{align*} Puisque $x$ et $y$ sont des entiers naturels, on obtient le couple $(1,2)$.
  2. Deuxième méthode :
    Le couple $(0,0)$ est une solution évidente de l'équation. Cherchons les autres solutions en supposant $(x,y) \neq (0,0)$.
    En posant $d=x\land y$, on a $d \neq 0$, et il existe deux entiers $a$ et $b$ tels que : \[ \begin{cases} x=da \\ y=db \\ a\land b=1 \end{cases} \]
    1. En substituant dans l'équation initiale, on trouve : \begin{align*} (da)^2((da)^2+7) &= db(2da+db) \\ d^2a^2(d^2a^2+7) &= d^2b(2a+b) \end{align*} Puisque $d \neq 0$, on peut diviser les deux membres par $d^2$ pour obtenir : \[ a^2(d^2a^2+7) = b(2a+b) \qquad (3) \]

    2. L'équation $(3)$ montre que $a$ divise le produit $b(2a+b)$. Comme $a\land b=1$, on déduit (par le lemme de Gauss) que : \[ a \mid (2a+b) \] Ce qui équivaut à dire que la fraction suivante est un entier : \[ \dfrac{2a+b}{a} = 2+\dfrac{b}{a} \] Cela implique que $a$ divise $b$. Or, comme $a$ et $b$ sont premiers entre eux, on a nécessairement : \[ a=1 \]

    3. L'équation $(3)$ s'écrit alors en remplaçant $a$ par $1$ : \begin{align*} d^2+7 &= b(2+b) \\ d^2+7 &= 2b+b^2 \\ d^2+8 &= b^2+2b+1 \\ d^2+8 &= (b+1)^2 \qquad (4) \end{align*}

    4. L'équation $(4)$ peut être factorisée de la manière suivante : \[ (b+1)^2-d^2 = 8 \implies ((b+1)-d)((b+1)+d) = 8 \] En remarquant que :
      • $((b+1)-d) < ((b+1)+d)$ (car $d>0$)
      • $((b+1)-d)$ et $((b+1)+d)$ ont la même parité.
      • Leur produit est pair (égal à $8$), donc les deux facteurs sont pairs.
      Le seul système possible est : \begin{align*} \begin{cases} (b+1)-d &= 2 \\ (b+1)+d &= 4 \end{cases} &\implies \begin{cases} 2(b+1) &= 6 \\ 2d &= 2 \end{cases} \\ &\implies \begin{cases} b+1 &= 3 \\ d &= 1 \end{cases} \\ &\implies \begin{cases} b &= 2 \\ d &= 1 \end{cases} \end{align*} On retrouve alors les valeurs de $x$ et $y$ : \[ \begin{cases} x = da = 1 \times 1 = 1 \\ y = db = 1 \times 2 = 2 \end{cases} \] Conclusion :
      L'ensemble des solutions de l'équation $(E)$ dans $\mathbb{N}^2$ est : \[ S = \{(0,0), (1,2)\} \]