First of all, we note that \(x\) and \(y\) cannot be both squares. Indeed, suppose that \(x\) is a square. Since \(xy+x=x(y+1) \) is a square, \(y+1\) must also be a square. But then (since \(y >0\) ) \(y \) cannot be a square.
Suppose now that \(x\) and \(y\) are positive integers such that both \( x(y+1) \) and \( y(x+1) \) are squares. We can write \( x = au^2\) and \(y= bv^2\) for some positive integers \(a, b, u, v\), where \(a\) and \(b\) are squarefree. We need to show that \(a =1 \) or \(b =1 \).
Since \(au^2(bv^2 +1) \) and \(bv^2(au^2 +1) \) are perfect squares, there exist positive integers \(s \) and \(t \) such that \(a\left( {b{v^2} + 1} \right) = {\left( {as} \right)^2}\) and \(b\left( {a{u^2} + 1} \right) = {\left( {bt} \right)^2}\). It follows that \(\boxed{a{s^2}  b{v^2} = 1}\) and \(\boxed{a{u^2}  b{t^2} =  1}\). Note that \(\gcd \left( {a,b} \right) = 1\).
Let \(z = st  uv\). Then, we have:
\(\left( {b{v^2}} \right)\left( {a{u^2}} \right) = ab{\left( {st  z} \right)^2} \Leftrightarrow \left( {a{s^2}  1} \right)\left( {b{t^2}  1} \right) = ab{\left( {st  z} \right)^2} \Leftrightarrow \)
\( \Leftrightarrow ab{s^2}{t^2}  a{s^2}  b{t^2} + 1 = ab{s^2}{t^2}  2abstz + ab{z^2} \Leftrightarrow \)
\( \Leftrightarrow a{s^2} + b{t^2} + ab{z^2}  1 = 2abstz\) \( \bf \color{red} {(1)} \).
If we set \(m = as^2 \), \( n = bt^2 \) and \(p = abz^2 \), then \( \bf \color{red} {(1)} \) is equivalent to \[\boxed{m + n + p  1 = 2\sqrt {mnp} } \bf \color{red} {(\bigstar)}.\] Claim : At least one of the numbers \(m, n \) and \( p\) is a perfect square.
It follows from the Claim that at least one of the numbers \(a, b\) and \(ab\) is a square. But \(a\) and \(b\) are squarefree and \(\gcd \left( {a,b} \right) = 1\), so \(a =1 \) or \( b = 1 \), answering the first part of the problem.
It remains to prove the Claim. Rewrite equation \(\bf \color{red} {(\bigstar)}\) as \({\left( {m + n + p  1} \right)^2} = 4mnp \bf \color {red} {(2)}\).
Let \(\alpha = 2m  1\), \( \beta = 2n 1 \) and \( \gamma = 2p1\). Then, equation \(\bf \color {red} {(2)}\) is equivalent to
\({\left( {\alpha + \beta + \gamma + 1} \right)^2} = 2\left( {\alpha + 1} \right)\left( {\beta + 1} \right)\left( {\gamma + 1} \right) \Leftrightarrow \)
\( \Leftrightarrow {\alpha ^2} + {\beta ^2} + {\gamma ^2}  2\alpha \beta \gamma = 1 \bf \color {red} {(3)}\).
Now equation \(\bf \color {red} {(3)}\) is equivalent to each of the equations
\(\left( {{\alpha ^2}  1} \right)\left( {{\beta ^2}  1} \right) = {\left( {\alpha \beta  \gamma } \right)^2}\)
and
\(\left( {{\beta ^2}  1} \right)\left( {{\gamma ^2}  1} \right) = {\left( {\beta \gamma  \alpha } \right)^2}\),
from which we infer that there exist nonnegative integers \(d, e, f, g\), with \( d \) not a perfect square, such that
\({\alpha ^2}  1 = d{e^2}\), \({\beta ^2}  1 = d{f^2}\), \({\gamma ^2}  1 = d{g^2}\),
\(\left {\alpha \beta  \gamma } \right = def\) and \(\left {\beta \gamma  \alpha } \right = dfg\).
Now we invoke the theory of the Pell equation. Let \(\left( {{x_1},{y_1}} \right)\) be the fundamental solution of the Pell equation \({X^2}  d{Y^2} = 1\). If we set \({z_1} = {x_1} + {y_1}\sqrt d \), then it is known that all its solutions \(\left( {{x_k},{y_k}} \right)\) (where \( k \ge 0\) ) are given by \[{x_k} = \frac{1}{2}\left( {z_1^k + z_1^{  k}} \right), \mbox{ } {y_k} = \frac{1}{{2\sqrt d }}\left( {z_1^k  z_1^{  k}} \right).\] It follows that there exist nonnegative integers \(k_1, k_2\) and \(k_3\), such that \[\alpha = \frac{1}{2}\left( {z_1^{{k_1}} + z_1^{  {k_1}}} \right),\beta = \frac{1}{2}\left( {z_1^{{k_2}} + z_1^{  {k_2}}} \right),\gamma = \frac{1}{2}\left( {z_1^{{k_3}} + z_1^{  {k_3}}} \right).\] Assume, without loss of generality, that \(m \ge n \ge p\). Then, we have \(k_1 \ge k_2 \ge k_3 \) and \(\alpha \beta  \gamma = def\). Hence, \[\gamma = \alpha \beta  def =\] \[ =\frac{1}{4}\left( {z_1^{{k_1}} + z_1^{  {k_1}}} \right)\left( {z_1^{{k_2}} + z_1^{  {k_2}}} \right)  \frac{1}{4}\left( {z_1^{{k_1}}  z_1^{  {k_1}}} \right)\left( {z_1^{{k_2}}  z_1^{  {k_2}}} \right) = \] \[ = \frac{1}{2}\left( {z_1^{{k_1}  {k_2}} + z_1^{  \left( {{k_1}  {k_2}} \right)}} \right).\] We conclude that \(k_3 = k_1  k_2\), so not all of the numbers \(k_1, k_2\) and \(k_3\) are odd. Suppose that \(k_1\) is even. Then, we can write \[m = \frac{{\alpha + 1}}{2} = {\left[ {\frac{1}{2}\left( {z_1^{\frac{{{k_1}}}{2}} + z_1^{  \frac{{{k_1}}}{2}}} \right)} \right]^2},\] so \(m\) is a perfect square and the Claim follows.
To answer the second question. i.e. to characterize all pairs \( (x,y) \) of positive integers with the property that both \(x(y+1) \) and \( y(x+1) \) are squares, we may proceed as follows: Without loss of generality, we may assume that \( x \) is a square. Then, \(y+1 \) will also be a square, so there exist positive integers \( r\) and \(s\) such that \(x = r^2 \) and \(y+1 = s^2\). We also have that \(y(x+1) = (r^2 +1)(s^2 1) \) is a square. Let \(d = \gcd \left( {{r^2} + 1,{s^2}  1} \right)\). Then, there exist positive integers \(t\) and \(u\) such that
\({r^2} + 1 = d{t^2} \Longleftrightarrow r^2  dt^2 =1\)
and
\({s^2}  1 = d{u^2} \Longleftrightarrow s^2du^2 =1\).
Thus, the pair \((r,t) \) of positive integers is a solution of the negative Pell equation \(X^2  dY^2 = 1 \). Depending on \(d\), this has either no solutions or infinitely many solutions.
According to the theory of the Pell equation, if \(d \) is such that the above negative Pell equation is solvable in the positive integers, then we can find positive integers \(r,s,t \) and \(u\) such that \({r^2} + 1 = d{t^2}\) and \({s^2}  1 = d{u^2} \). Then the pair \[\left( {x,y} \right) = \left( {{r^2},d{u^2}} \right)\] has the desired property.
One could possibly characterize (via continued fractions) all \(d\) such that the negative Pell equation \(X^2  dY^2 = 1 \) has (infinitely many) solutions and therefore obtain a more precise answer to the second part of the problem.
