Let $x$ and $y$ be positive integers such that $xy \mid x^2+y^2+1$.

Use the very tricky technique called Vieta Jumping.
The idea is considering a polynomial $f(x,y)$ that is quadratic in both $x$ and $y$, with integer coefficients and symmetrical (that is $f(x,y)=f(y,x)$. We have that if $f(x,y)$ has some property when $x,y$ are integers and we want to prove something regarding $x$ and $y$. Suppose that some pair $x_1,y_1$ of integers satisfies the property, since $f$ is symmetrical, we can suppose WLOG that $x_1>y_1$ (the case $x_1=y_1$ is usually easy).

Recall the vieta formulas:
If $z_1$ and $z_2$ are the roots of $x^2+bx+c=0$, then $z_1z_2=c$ and $z_1+z_2=-b$.
Those formulas are very useful, particularly the last one, since it is a simple sum.

Now since $f(x,y)$ is quadratic in $x$, we apply the vieta formulas in $x_1$ and we find some integer $x_0$ with $x_0<y_1$ that satisfy the same property, Now we do the same with $y_1$ and find another integer $y_0$ with $y_0<x_0$ that also satisfy the property. Continuing this way we get a pair $(a,b)$ of integers that satisfy the property with $a$ and $b$ really small (like $a=1$). It's easy to prove what we want when the integers are small. Now since all these pairs were satisfying the same property, what we proved about $(a,b)$, also applies to the initial $(x_1,y_1)$.

Well, that was kinda long. I hope i have explained the main point. Try to use this on the problem and then back and post your results :)

Suppose $xy\mid x^2+y^2+1$ and let $t=\displaystyle\frac{x^2+y^2+1}{xy}$ such that $t\in\mathbb{N}$.

Construct the set, $$S=\left\{(x,y)\in\mathbb{N}\times\mathbb{N} : \frac{x^2+y^2+1}{xy}=t\in\mathbb{N}\right\}$$

We deduce that $\displaystyle\frac{x^2+y^2+1}{xy} \ge 3$ because $\displaystyle\frac{x^2+y^2+1}{xy}<3$ implies $x^2+y^2+1 \le 2xy \le x^2+y^2$ which is clearly a contradiction. Now fix $t$ and suppose that $t>3$. Since $S\neq \varnothing$, we can choose $(a,b)\in S$ such that $a+b$ is minimal and $t>3$. WLOG assume $a\ge b>0$. Let us consider the quadratic

$$p(w)=w^2-tbw+b^2+1=0$$

It follows that $a$ is a solution since $(a,b)\in S$ and hence satisfies the quadratic equation, that is $a$ is a root. By applying Vieta's formulas we obtain the other root $c$. Hence $a+c=tb$ and $ac=b^2+1$. Since $c=tb-a$ we have $c\in\mathbb{Z}$. Now it remains to prove that $(a,c)\in S$.

To this end suppose $c<0$. It follows that $$\displaystyle 0<a^2+ac+1-3ab=a^2+ac-\frac{3c}{t}(a+c)<0$$ which is clearly a contradiction. It also immediately follows that $c\neq 0$ as this implies $b^2+1=0$. Therefore $c\in\mathbb{N}$ and $(a,c)\in S$.

Now we show that this $c$ contradicts the minimality of $a$, that is $c<a$. Suppose $c>a$ so it follows that $a+1\le c$. But from Vieta's equations we obtain $\displaystyle a+1\le c=\frac{b^2+1}{a}\le a+\frac{1}{a}$ which is impossible since this inequality holds in $\mathbb{N}$ if and only if $a=1$ and hence $a=b=1$ implying $t=3$ which contradicts our assumption that $t>3$. Therefore $c\le a$. But if $c=a$ then this implies that $\displaystyle a^2=b^2+1>\frac{9}{4}b^2$ which again is a contradiction. Hence we conclude that $c<a$ and as a result $c+b$ contradicts the minimality of $a+b$. Hence $t=3$

I found this on Wolfram under the Fibonacci Number page:

http://mathworld.wolfram.com/FibonacciNumber.html

Catalan's Identity:

$F(n)^2 - F(n-r)F(n+r) = (-1)^{(n-r)}F(r)^2$

For $r = 1$, n is even

$F(n)^2 - F(n-1)F(n+1) = -1$

Replace F(n) using $F(n) = F(n+1) - F(n-1)$ and you get

$(F(n+1) - F(n-1))^2 - F(n-1)F(n+1) = -1$

or

$F(n+1)^2 + F(n-1)^2 + 1 = 3F(n-1)F(n+1)$

is of the form:

$x^2 + y^2 + 1 = 3xy$

I do not believe it shows that only Fibonacci numbers are solutions. The relationship that I was looking for on my other solution uses Catalan's Identity, with r = 2.