Finding $(a,b)\in\mathbb{N}^2$ such that $\dfrac{a^2+b^2+1}{a+b} \in \mathbb{N}$. A pair $(a,b)\in\mathbb{N}^2$ is called good if $a < b$ and
$$\frac{a^2+b^2+1}{a+b}\in\mathbb{N}.$$
I think I've shown that there are infinitely many good pairs. However, the family of good pairs that I've found is not all of them, and the argument is not as elementary as I'd like.
Observation: we can check small values to show that good pairs do indeed exist, e.g. $(1,2),\, (4,7),\,(5,12),\, (11,16)$.
To construct an infinite family of good pairs, suppose we have
$$\frac{a^2+b^2+1}{a+b}=k\in\mathbb{N},$$
then $a^2+b^2+1 - k(a+b)=0$, and completing the square for $a$ and $b$ gives the equation for a circle:
$$\left(a-\frac{k}{2}\right)^2+\left(b-\frac{k}{2}\right)^2 = \frac{k^2}{2}-1$$
$$\implies a = \frac{k}{2} + \sqrt{\frac{k^2}{2}-1-\left(b-\frac{k}{2}\right)^2}.$$
Let $k=2j$ for some $j\in\mathbb{N}$ to eliminate the fractions, obtaining $a=j+\sqrt{j^2-1-b^2+2jb}$. In order to have $a\in\mathbb{N}$, we must have $j^2-1-b^2+2bj=n^2$ for some $n\in\mathbb{N}\cup\{0\}$. We can view this as a quadratic in $b$:
$$b^2-2jb+(n^2+1-j^2)=0.$$
Since we require $b\in\mathbb{N}$, the discriminant $\Delta$ must be of the form $\Delta = 4m^2$ for some $m\in\mathbb{N}\cup\{0\}$. Thus we have the following expressions for $a$ and $b$:
$$b=\dfrac{2j\pm\sqrt{\Delta}}{2} = j+m, \qquad a = j + \sqrt{n^2} = j+n.$$
Considering the discriminant:
$$\Delta = (-2j)^2 - 4(n^2+1-j^2) = 4m^2$$
$$\implies 2j^2-1 = m^2+n^2.$$
If $a<b$ we must have $n<m$. In particular, we can assume $n=0$ without restricting $m$. Then we have $m^2-2j^2=-1$, and so we have a correspondence between certain good pairs $(a,b)$ and elements of norm $-1$ in $\mathbb{Q}(\sqrt{2})$. For $x+y\sqrt{2}\in\mathbb{Z}[\sqrt{2}]$, define $N(x+y\sqrt{2})=x^2-2y^2$, and observe that this function is multiplicative. We have $N(1\pm\sqrt{2})=-1$, and thus $N((1\pm\sqrt{2})^{2r-1})=(-1)^{2r-1}=-1$ for any $r\in\mathbb{N}$. Let $m+j\sqrt{2} = (1+\sqrt{2})^{2r-1}$ and so $m-j\sqrt{2} = (1-\sqrt{2})^{2r-1}$. We can combine these expressions to obtain the following general formulas for values of $m$ and $j$:
$$m_r = \frac{(1+\sqrt{2})^{2r-1}+(1-\sqrt{2})^{2r-1}}{2},\qquad j_r = \frac{(1+\sqrt{2})^{2r-1}-(1-\sqrt{2})^{2r-1}}{2\sqrt{2}}.$$
Then we have $a=j+n=j$ and $b=m+j$, thus we have an infinite family of good pairs $(a_r,b_r)$ for $r\in\mathbb{N}$ given by:
$$a_r = \frac{(1+\sqrt{2})^{2r-1}-(1-\sqrt{2})^{2r-1}}{2\sqrt{2}},\qquad b_r = \frac{(1+\sqrt{2})^{2r}-(1-\sqrt{2})^{2r}}{2\sqrt{2}}.$$
The first few values of $r$ gives the good pairs $(1,2),\,(5,12),\,(29,70),\,(169,408),(985,2378)$ with corresponding values of $k=2,10,58,338, 1970$.
I seek advice both on this proof itself (i.e. does it work?) and also in seeking a more elementary argument (i.e. without requiring knowledge about $\mathbb{Q}(\sqrt{2})$ and the like) - it doesn't necessarily need to be constructive, in fact it would be interesting to know if there exists a more concise argument to that effect.
Regarding finding good pairs themselves, comparing the list constructed in my solution to thta of the observation, it's clear that this solution 'misses' some pairs - and this comes a no surprise since we disregard the cases that $k$ is odd, and $n>0$.
If $k=2j-1$, then we can use a similar method to show $j^2-j-b^2+2jb-b=n^2-n$ for some $n\in\mathbb{N}$, but I haven't been able to make any progress from this. The case for $k$ even and $n>0$ is essentially asking "when can $2j^2+1$ be written as the sum of two squares", which I'm not sure quite how to answer - using known results about sums of two square seems difficult as $2j^2+1$ has no general factorisation.
Are there totally different methods that can be used to construct more families of good pairs? And what about any other 'structure' to good pairs? It seems interesting that units in $\mathbb{Z}[\sqrt{2}]$ give rise to the solutions, can a more general statement be made?

$a^2-b^2$ = $(a-b)(a+b)$ is divisible by a+b, so $a^2+b^2+1$ is divisible by a+b if and only if $2a^2+1$ is.
To find all solutions with a, b >= 1:
Iterate through a = 0, 1, 2, …
Let z = 2a^2+1
Determine all divisors d of z. 
For each d >= a, (a, d-a) is a solution. 

Since z is divisible by z, we have the obvious b = z - a for each a. For example, a = 5, z = 51, b = 46, and $5^2 + 46^2 + 1 = 2142 = 51 \cdot 42$. Other divisors of z are 1, 3 and 17, leading to be = 17-5 = 12, and $5^2 + 12^2 + 1 = 170 is divisible by a+b = 17.
The values z have on average fewer divisors than random integer of the same size, since they are never divisible by 2, 5, 7, 11, or 13.
