<s>[INST] What are the non-negative integers of this form? What are $( \frac{x^2-1-my}{x}, \frac{y^2-1-nx}{y} ) \in  \mathbb{Z}^2_{\ge 0}$ for $(m,n) \in \mathbb{Z}^2_{\ge 0}$, $(x,y) \in \mathbb{R}^2_{\ge 1}$ and $mx+ny=xy$?  

Some investigations:   


*

*According to this answer of Ivan Neretin, assuming $x,y$ integers, the only possibility is $(0,1)$ with $(m,n) = (0,1)$ and $(x,y) = (1,2)$, up to permutation.

*If $x=y$, then $(m+n)x = x^2$, so $x=m+n$ positive integer.


In general $y=mx/(x-n)$ with $x>n$. Thus we need to consider
$$ ( \frac{x^2-1-m\frac{mx}{x-n}}{x}, \frac{(mx/(x-n))^2-1-nx}{mx/(x-n)}  ) $$ $$ = (\frac{(x^2-1)(x-n)-m^2x}{x(x-n)},\frac{m^2x^2-(x-n)^2(1+nx)}{mx(x-n)}) \in  \mathbb{Z}^2_{\ge 0} $$ 
Let $\alpha=\alpha_{n,m}$ be the root of $$P_{n,m}(x):=(x^2-1)(x-n)-m^2x$$ when $x>n$ (show unicity). Then $$m^2\alpha^2-(\alpha-n)^2(1+n\alpha) = (\alpha^2-1)(\alpha-n)\alpha -(\alpha-n)^2(1+n\alpha) $$ $$ = (\alpha-n)((\alpha^2-1)(\alpha-n)-\alpha+n^2\alpha) = (m^2+n^2-1)\alpha(\alpha-n)$$ 
Let $\beta = \beta_{n,m}$ be the root of $$Q_{n,m}(x):=m^2x^2-(x-n)^2(1+nx)$$ when $x>n$ (show unicity). Note that if $n=m$ then $\beta = n \alpha/(\alpha-n)$.
We can assume (by symmetry) that $\alpha \le \beta$.
Let $(a,b) = ( \frac{x^2-1-my}{x}, \frac{y^2-1-nx}{y} )$. Then for $x=\alpha$ we get $$(a,b) = (0,\frac{m^2+n^2-1}{m}),$$ and by symmetry, for $x=\beta$, we get $$(a,b) = (\frac{m^2+n^2-1}{n},0).$$
The problem reduces to find all the $x \in [\alpha, \beta]$ such that $P_{n,m}(x)$ an $Q_{n,m}(x)$ are integers.    
Remark: If $n |(m^2-1)$ or $m |(n^2-1)$, then above is at least one solution.
Example: $m=n=1$, then $\alpha=1.80193773580484...$ and $\beta=2.24697960371747...$.
Let us plot $P_{1,1}$ and $Q_{1,1}$:

It follows that $(0,1)$ and $(1,0)$ are the only solutions if $m=n=1$.  
What are all the others solutions in general?

An other approach 
Take $(p,q) = ( \frac{x^2-1-my}{x}, \frac{y^2-1-nx}{y} )$. Then: 
$$\left\{
\begin{array}{ll}
x^2 - px-1-my = 0 \\
y^2 -qy-1-nx = 0
\end{array}
\right.$$
Because $(x,y) \in \mathbb{R}^2_{>1}$, it follows that:   
$$\left\{
\begin{array}{ll}
2x = p + (p^2+4(1+my))^{1/2}\\
2y = q + (q^2+4(1+nx))^{1/2}
\end{array}
\right.$$
So $$ 2x = p + [p^2+4+2m(q + [q^2+4(1+nx)]^{1/2})]^{1/2} $$
Then
$$ [(2x-p)^2-p^2-4-2mq]^2 = q^2+4(1+nx) $$
Thus
$$ (4x^2-4xp-4-2mq)^2 = q^2+4(1+nx) $$ 
We got a polynomial of degree $4$, allowing us to get at most $4$ roots $x_i(m,n,p,q)$ with $i \le 4$. Then $$2y_i(m,n,p,q) = q + (q^2+4(1+nx_i(m,n,p,q)))^{1/2}$$
It remains to solve the Diophantine equation $$mx_i(m,n,p,q)+ny_i(m,n,p,q)=x_i(m,n,p,q)y_i(m,n,p,q).$$ 
 [/INST] A pair $(p,q) = ( \frac{x^2-1-my}{x}, \frac{y^2-1-nx}{y} )$ is achievable if and only if $p^2+q^2+4$ is the sum of two squares.

Proof. In one direction, we assume that the following three polynomials have a common zero: 
\begin{split}
f(x,y)&:=x^2-1-my-px,\\
g(x,y)&:=y^2-1-nx-qy, \\ 
h(x,y)&:=xy-mx-ny.
\end{split} Then eliminating $x,y$ via computing resultants, we have with necessity:
$$0 = \mathrm{Res}_y(\mathrm{Res}_x(f,g),\mathrm{Res}_x(g,h)) = -(m^2-mq+n^2-np-1)^3 n^4,$$
implying that $m^2-mq+n^2-np-1=0$, which is equivalent to 
$$(\star)\qquad p^2 + q^2 + 4 = (2n-p)^2 + (2m-q)^2.$$ 
That is, $p^2+q^2+4$ is the sum of two squares.
In the reverse direction, we assume that $p^2 + q^2 + 4 = a^2 + b^2$ for some nonnegative integers $a,b$. Without loss of generality we have $a\equiv p\pmod{2}$ and $b\equiv q\pmod{2}$. We set $n:=\frac{a+p}2$ and $m:=\frac{b+q}2$, which makes them satisfy $(\star)$. 
Then we set $x_0>n$  be a zero of the polynomial:
$$F(x):=x^3-(p+n)x^2+(pn-1-m^2)x+n.$$
Such zero does exist since $F(n) = -m^2n<0$ while $F(x)\to+\infty$ as $x$ grow. Then, we set $y_0:=\tfrac{mx_0}{x_0-n}$ to get $h(x_0,y_0)=0$. Since $F(x) = (x-n)f(x,y)+mh(x,y)$, we also have $f(x_0,y_0)=0$, which together with $(\star)$ further implies that $g(x_0,y_0)=0$. That is, for $(x,y)=(x_0,y_0)$ we get $( \tfrac{x^2-1-my}{x}, \tfrac{y^2-1-nx}{y} ) = (p,q)$. QED
P.S. Similarly to $x_0$, we also have that $y_0$ is a zero of the polynomial:
$$G(y):= y^3-(q+m)y^2+(qm-1-n^2)y+m.$$
which satisfies $G(y)=(y-m)g(x,y) + nh(x,y)$.
</s>