On the diophantine equation $x^{m-1}(x+1)=y^{n-1}(y+1)$ with $x>y$, over integers greater or equal than two

I don't know if the following diophantine equation (problem) is in the literature. We consider the diophantine equation $$x^{m-1}(x+1)=y^{n-1}(y+1)\tag{1}$$ over integers $x\geq 2$ and $y\geq 2$ with $x>y$, and over integers $m\geq 2$ and $n\geq 2$. These are four integral variables $x,y,m$ and $n$. The solutions that I know for the problem $(1)$ are two, the solution $(x,y;m,n)=(3,2;2,3)$ and $(98,21;2,3)$.

Question 1. Do you know if this problem is in the literature? Alternatively, if this problem isn't in the literature can you find more solutions?

If the equation or problem $(1)$ is in the literature please refer it answering this question as a reference request, and I try to search and read the statements for new solutions from the literature. In other case compute more solutions or add upto what uppers limits you got evidence that there aren't more solutions.

Question 2. I would like to know what work can be done with the purpose to know if the problem $(1)$ have finitely many solutions $(x,y;m,n)$. I mean what relevant reasonings or heuristics you can to deduce with the purpose to study if the problem have finitely solutions.

If this second question is in the literature, please refer the literature answering this question as a reference request, and I try to search and read the statements from the literature.

This answer is open for improvement. Feel free to use its results answering the question.

Conjecture 1. There are no solutions when $m|n$.

For each integer $p\ge 2$ and each real $t$ put $f_p(t)=t^{p-1}(t+1)$. In order to prove the conjecture for each integer $y\ge 2$ we hope to find $x’$ and $x’’$ such that $f_m(x’)<f_n(y)<f_m(x’’)$ but there are no integer values between $x’$ and $x’’$.

Claim 2. There are no solutions when $n=2m$.

Proof. We claim that $x’=y^2+\tfrac ym-\tfrac 2m$ and $x’’=x’+\tfrac 1m$ fit. Indeed, $$f_n(x’’)=\left(y^2+\frac ym-\frac 1m\right)^{m}+\left(y^2+\frac ym-\frac 1m\right)^{m-1}>$$ $$y^{2m}+{m\choose 1} y^{2m-2}\left(\frac ym-\frac 1m\right)+ y^{2m-2}= y^{2m}+ y^{2m-1}=f_n(y).$$

I have a draft proof that $x’$ fits based on Bernoulli’s inequality, but it is cumbersome.

Proposition 3. There are no solutions when $m=2$ and $n=6$.

We claim that $x’=y^3+\tfrac {y^2}2-\tfrac y8-\tfrac 12$ and $x’’=x’+\tfrac 18$ fit. Indeed, we can check (by Mathcad) that

$$f_m(x’’)-f_n(y)=\frac 1{2^6}(y-1)\left(8y^2+17y+15\right)>0$$ and $$f_n(y)- f_m(x’)=\frac 1{2^6}\left(8y^3-y^2+16\right)>0.$$

Proposition 4. There are no solutions when $m=3$ and $n=9$.

We claim that $x’=y^3+\tfrac {y^2}3-\tfrac {y}{9}-\tfrac 1{3}$ and $x’’=x’+\tfrac 1{9}$ fit. Indeed, we can check (by Mathcad) that

$$f_m(x’’)-f_n(y)=\frac 1{3^6}(y-1)\left(108y^5+270y^4+252y^3+17y^2-52y-28\right)>0$$ and

$$f_n(y)- f_m(x’)= \frac 1{3^6} \left(135y^6-9y^4+244y^3+81y^2-27y-54\right)>0.$$

Proposition 5. There are no solutions when $m=4$ and $n=12$.

We claim that $x’=y^3+\tfrac {y^2}4-\tfrac {3y}{32}-\tfrac 7{32}$ and $x’’=x’+\tfrac 1{32}$ fit. Indeed, we can check (by Mathcad) that

$$f_m(x’’)-f_n(y)=$$ $$\frac 1{2^{20}}\left(32768y^9+182272y^8-3072y^7-401408y^6-182880y^5+47889y^4+100008y^3+17496y^2-7776y-5616\right)>0$$ and

$$f_n(y)- f_m(x’)=$$ $$\frac 1{2^{20}}\left(98304y^9-83968y^8-9216y^7+403456y^6+193248y^5-49329y^4-124660y^3-23254y^2+9996y+8575\right)>0.$$

We see that for bigger $m$ and $n$, expressions for $x’$ and $x’’$ get more and more complicated. So we have

Working problem 6. Is there a pattern for $x’$ and $x’’$ and, if so, then can we prove the conjecture from this pattern?

My approach is to restrict the exponents to $m = 2, n = 3$, because both of your solutions have these values, and due to $m - 1 = 1$, this could makes the question different, as the LHS becomes simply a product of two consecutive numbers.
This approach of mine might not be very useful, as I thought I can come up with something that proves that there are infinitely many solutions for this type of equality, but at the end you see it might be the other way round, your two solutions might be the only ones.
I will use the identity $a^b - 1 = (a - 1)(a^{b - 1} + a^{b - 2} + ... + a^2 + a + 1)$ in my approach. $\tag{*}$
We have $x(x + 1) = y^2(y + 1) = y^3 + y^2$
and to be able to use $(*)$, subtract $2$ from both sides:
$x(x + 1) - 2 = y^3 + y^2 - 2$
$x^2 + x - 2 = y^3 - 1 + y^2 - 1 = (y - 1)(y^2 + y + 1) + (y - 1)(y + 1)$
Now I would like to have an expression on both sides so that I can piecewise make them equal:
$x^2 + x - 2 = (y - 1)(y^2 + 2y + 1) + y - 1$
$x^2 + x - 2 = (y - 1)(y + 1)^2 + (y + 1) - 2$
$x^2 + x = (y - 1)(y + 1)^2 + (y + 1)$
Now we had on both sides pairwise equal terms of $x = (y + 1)$, if $(y - 1) = 1$. So one solution is $y = 2, x = (y + 1) = 3$. This is what you got. But this is not giving your other solution. Something is needed to make it more general.
My idea is that as we already have an expression on the RHS that has the right exponents, we need to get rid of the coefficient $(y - 1)$. Therefore we change it to
$x^2 + x = (y - 1)(y + 1)^2 + (y + 1) = (y + 1)^2 + (y + 1) + (y - 2)(y + 1)^2$
and so we must have $(y - 2)(y + 1)^2 = 0$ to obtain $x = (y + 1)$. But this only give $y = -1$ besides $y = 2$, so nothing useful.
So let us change the delta between $x$ and $y$, and instead of having $x = (y + 1)$, let us change to $x = (y + \alpha)$.
We need the following:
$(y + \alpha) = (y + 1) + (\alpha - 1) \implies (y + 1) = (y + \alpha) - (\alpha - 1)$
$(y + \alpha)^2 = y^2 + 2\alpha y + \alpha^2 = y^2 + 2y + 1 + (2\alpha - 2)y + (\alpha^2 - 1) = (y + 1)^2 + (2\alpha - 2)y + (\alpha^2 - 1) \implies (y + 1)^2 = (y + \alpha)^2 - (2\alpha y - 2y + \alpha^2 - 1)$
Now we have
$x^2 + x = (y + 1)^2 + (y + 1) + (y - 2)(y + 1)^2 = (y + \alpha)^2 + (y + \alpha) - (2\alpha y - 2y + \alpha^2 - 1) - (\alpha - 1) + (y - 2)((y + \alpha)^2 - (2\alpha y - 2y + \alpha^2 - 1))$
and so we must have $ - (2\alpha y - 2y + \alpha^2 - 1) - (\alpha - 1) + (y - 2)((y + \alpha)^2 - (2\alpha y - 2y + \alpha^2 - 1)) = 0$ to obtain $x = (y + \alpha)$.
After a lot of transformations one gets a quadratic for $\alpha$:
$\alpha_{1, 2} = \alpha^2 + (2y + 1)\alpha + y - y^3$
from which $\alpha$ based on the quadratic formula is
$\alpha_{1, 2} = \frac{-(2y + 1) \pm \sqrt{4y^3 + 4y^2 + 1}}{2}$
I thought I can come up with something useful for the expression under the square root, but I could not determine for what values of $y$ it gives a square number. When it does, $\alpha$ is an integer, as the number under the square root is odd, therefore its square root is odd, and $(2y + 1)$ is odd, so the $\pm$ of these are even. For $y = 2$ it gives $\alpha = 1$ so $x = 3$, which is your first solution and my previous one, and for $y = 21$ it gives $\alpha = 77$ so $x = 98$ which is your other solution. I brute-forced the expression under the square root up to one million, but there are no other solutions.
Of course even if it is provable that there are no other solutions for this expression to be a square, other solutions might exist for the original equality, as it only related to the piecewise equality logic applied.

In this answer we address

Question 2. I would like to know what work can be done with the purpose to know if the problem $(1)$ have finitely many solutions $(x,y;m,n)$. I mean what relevant reasonings or heuristics you can to deduce with the purpose to study if the problem have finitely solutions.

Let $m$ and $n$ are fixed. According to [G], Diophantine equations with two variables of degree greater than two have infinitely many (integer) solutions only in very rare cases. In particular, by a special and very complicated method K. Zigel’ (Siegel?) showed the following

Theorem. Let $P(x,y)$ be a irreducible polynomial of two variables with integer coefficients of a total degree greater than two (that is, $P(x,y)$ contains a monomial $ax^ky^s$, where $k+s>2$). (The irreducibility of $P(x,y)$ means that it cannot be represented as a product of two non-constant polynomials with integer coefficients). If an equation $P(x,y)=0$ has infinitely many integer solutions $(x,y)$ then there exist an integer $r$ and integers $a_i$, $b_i$ for each $-r\le i\le r$ such that if in the equation $P(x,y)=0$ we make a substitution $x=\sum_{i=-r}^r a_it^i$ and $y=\sum_{i=-r}^r b_it^i$ then we obtain an identity.

For your question $P(x,y)=y^n+y^{n-1}-x^m-x^{m-1}$. In User2020201’s answer is shown that $m<n$. GreginGre showed that $P(x,y)$ is irreducible. It looks plausible and easy to prove that there is no above substitutions making $P(x,y)$ identity. But I am not a specialist in this branch of mathematics, so I asked a separate question for help.

References

[G] Gel’fand A.O. Solutions of equations in integer numbers, 3-rd edn., Moscow, Nauka, 1978, in Russian.