This is an additional remark, not really an answer to the question. I shall not provide a proof either, but the technique I used to get this result is basically Vieta Jumping. I have to post as an answer because the comment is too long.
It is well known (at least you can find a proof of this extraordinary claim in the link given above) that, if $a,b\in\mathbb{N}_0$ are such that $k:=\frac{a^2+b^2}{ab+1}$ is an integer, then $k=m^2$ for some $m\in\mathbb{N}_0$. For a fixed $m$, every solution $(a,b)\in\mathbb{N}_0\times\mathbb{N}_0$ to $\frac{a^2+b^2}{ab+1}=m^2$ is either of the form $\left(x_r(m),x_{r+1}(m)\right)$ or of the form $\left(x_{r+1}(m),x_r(m)\right)$, where $\left(x_r(m)\right)_{r=0}^\infty$ is the sequence defined by $x_0(m)=0$, $x_1(m)=m$, and $$x_r(m)=m^2\cdot x_{r-1}(m)-x_{r-2}(m)$$ for all integers $r\geq 2$. Note that
$$x_r(m)=\frac{m}{\sqrt{m^4-4}}\left(\left(\frac{m^2+\sqrt{m^4-4}}{2}\right)^r-\left(\frac{m^2-\sqrt{m^4-4}}{2}\right)^r\right)$$
for all $r\in\mathbb{N}_0$.
Note that, for $m\geq 2$, $\left(x_r(m)\right)_{r=0}^\infty$ is strictly increasing. Hence, if $r\geq 1$ and $x_r(m)\mid x_{r+1}(m)$, then by the relation $x_{r+1}(m)=m^2\cdot x_{r}(m)-x_{r-1}(m)$, we have that $x_{r}(m)\mid x_{r-1}(m)$. This can occur only when $r=1$ (where $x_{r-1}(m)=x_0(m)=0$). Thus, if $(a,b)=(d,n)$ is a solution in $\mathbb{N}\times\mathbb{N}$ with $d\mid n$ to $\frac{a^2+b^2}{ab+1}=m^2$ with $m\geq 2$, then $$(d,n)=\left(x_1(m),x_2(m)\right)=\left(m,m^3\right)\,.$$ For $m=1$, it is obvious that the only solution $(a,b)\in\mathbb{N}\times\mathbb{N}$ to $\frac{a^2+b^2}{ab+1}=m^2$ is $(a,b)=(1,1)$.