First I will show that $ a\ne b+1$. If instead $ a=b+1$, then we would have $ b+1=a|b^2+b+1$ $ b+1|b(b+1)=b^2+b$ $ b+1|(b^2+b+1)-(b^2+b)=1$ which is impossible since $ b\ge 1$ and hence $ b+1\ge 2$. Now, I will show the problem is equivalent to $ a^2 + a + b^2 + b + 1\equiv 0\mod ab$. We know that $ a$ and $ b$ must be relatively prime. Otherwise, $ \gcd(a,b)|a^2 + a$ so $ \gcd(a,b)\not| a^2 + a + 1$ if $ \gcd(a,b)>1$. Now $ a^2 + a + b^2 + b + 1\equiv 0\mod a$ and $ a^2 + a + b^2 + b + 1\equiv 0\mod b$, and combining these, since $ a$ and $ b$ are relatively prime, gives $ a^2 + b^2 + a + b + 1\equiv 0\mod ab$. This goes the other way to. If $ a^2 + b^2 + a + b + 1\equiv 0\mod ab$ then $ a^2+a+b^2+b+1\equiv 0\mod a\implies b^2+b+1\equiv 0\mod a$ and $ a^2+a+b^2+b+1\equiv 0\mod b\implies a^2 + a + 1\equiv 0\mod a$. So now we just need to solve $ a^2+a+b^2+b+1\equiv0\mod ab$. Let $ a^2 + a + b^2 + b + 1 = kab$ for some integer $ k$. Also, assume that $ a>b$ and that $ a,b>1$. Then $ a^2 - (kb-1)a + b^2 + b + 1 = 0$ This is a quadratic equation, so if $ a$ is a solution, then by Vieta's, $ a'=kb-1-a=\frac{b^2+b+1}{a}$ is also a solution. Since $ b^2+b+1$ and $ a$ are postive, $ a'$ must be positive. Now (by assumption that $ a>b$), $ a'=\frac{b^2+b+1}{a}<\frac{b^2+b+1}{b}=b+1+\frac{1}{b}$ $ a'\le b+1$. (because $ b>1$ by assumption, we know that $ 0<\frac{1}{b}<1$, but $ a'$ is an integer). Also, $ a'\ne b+1$ by above, so $ a'\le b$. So $ (a',b)$ is another solution if we already had $ (a,b)$. Since $ a>b$ and $ a'\le b$, our $ (a',b)$ solution is smaller. So, we know that we can always produce a smaller solution as long as $ a,b>1$ and $ a\ne b$. Now let's look at the cases where either $ a=b$ or one of $ a,b$ is $ 1$. $ a=b$. We know that $ a$ and $ b$ are relatively prime so for $ a$ to equal $ b$ we must have $ a=b=1$. $ a=1$ or $ b=1$. If $ a=1$ then $ b$ divides $ 1^2+1+1=3$, so the only pairs with $ 1$ in them are $ (1,3)$ and $ (1,1)$. So we have $ (1,1)$ and $ (1,3)$. Note than $ (1,3)$ can be reduced, as $ 5ab=1^2+1+3^2+3+1$, so $ k=5$ and $ a'=kb-1-a=5(1)-1-3=1$, so $ (1,3)$ is reduced to $ (1,1)$. Thus, all pairs besides $ (1,1)$ can be reduced in this manner. Since this reduction process does not change $ k$, and $ k=5$ for $ (1,1)$, we have $ k=5$ for all solutions. So since all solutions can be reduced to $ (1,1)$ by $ (a,b)\to (5b-a-1,b), (a>b)$ any solution can be obtained by the reverse procedure $ (a,b)\to (5b-a-1,b), (a\le b)$. Therefore if we define the sequence $ a_0=1,a_1=1$, and $ a_{n+2}=5a_{n+1}-a_{n}-1$, all solutions $ (a,b)$ are consecutive terms of this sequence.