Note that if $(m, n)$ is a solution then $(mn-1) \mid (n^3-1)$ so $(mn-1) \mid m^3(n^3-1)$. Also $(mn-1)\mid m^3n^3-1$ so $(mn-1) \mid m^3-1$, so $(n, m)$ is also a solution.
We may thus WLOG assume $m \geq n$. If $n=1$ then $\frac{n^3-1}{mn-1}=0$ so $(m, 1)$ is a solution for all positive integers $m>1$. I would reject $(1, 1)$ as that leads to $\frac{0}{0}$.
Otherwise $n \geq 2$. Note that $\frac{n^3-1}{mn-1} \equiv 1 \pmod{n}$ so we may write $\frac{n^3-1}{mn-1}=kn+1$, for some $k \in \mathbb{Z}$.
Now $m, n \in \mathbb{Z}^+$ so $0<\frac{n^3-1}{mn-1}=kn+1$. Since $n \geq 2$, this implies $k \geq 0$.
Since $m \geq n$, we have $kn+1=\frac{n^3-1}{mn-1} \leq \frac{n^3-1}{n^2-1}=\frac{n^2+n+1}{n+1}=n+\frac{1}{n+1}<n+1$. Thus $k<1$.
Therefore $k=0$, so $\frac{n^3-1}{mn-1}=1$, so $n^3-1=mn-1$, so (since $m, n>0$) $m=n^2$.
In conclusion, all positive integer solutions are given by $(m, 1)$ for $m>1$, $(1, n)$ for $n>1$, $(m, m^2)$ for $m>1$, and $(n^2, n)$ for $n>1$.
Here is an alternative approach. Let $x_n,y_n$, and $z_n$ be as in the argument given in the question; clearly $x_1=2$, and $y_1=z_1=1$. For $n\ge 1$ let $X_n,Y_n$, and $Z_n$ be the sets of $n$-digit numbers using only the digits $2,3,7$, and $9$ and congruent modulo $3$ to $0,1$, and $2$, respectively (so that $x_n=|X_n|$, $y_n=|Y_n|$, and $z_n=|Z_n|$). Finally, recall that an integer is congruent modulo $3$ to the sum of its digits.
Now let $n>1$, and let $k\in X_n\cup Y_n\cup Z_n$. Let $\ell$ be the $(n-1)$-digit number obtained by removing the last digit of $k$, and let $d$ be the last digit of $k$. If $d=3$ or $d=9$, then $k\equiv\ell\pmod3$; if $d=7$, then $k\equiv\ell+1\pmod3$; and if $d=2$, then $k\equiv\ell+2\pmod3$. Thus, $k\in X_n$ iff either $d\in\{3,9\}$ and $\ell\in X_{n-1}$, or $d=2$ and $\ell\in Y_{n-1}$, or $d=7$ and $\ell\in Z_{n-1}$. These three cases are exhaustive and mutually exclusive, so
$$x_n=|X_n|=2|X_{n-1}|+|Y_{n-1}|+|Z_{n-1}|=2x_{n-1}+y_{n-1}+z_{n-1}\;.$$
Similar reasoning shows that
$$y_n=2y_{n-1}+x_{n-1}+z_{n-1}$$
and
$$z_n=2z_{n-1}+x_{n-1}+y_{n-1}\;.$$
But clearly $x_{n-1}+y_{n-1}+z_{n-1}=4^{n-1}$, so these recurrences reduce to
$$\left\{\begin{align*}
x_n&=x_{n-1}+4^{n-1}\\
y_n&=y_{n-1}+4^{n-1}\\
z_n&=z_{n-1}+4^{n-1}\;.
\end{align*}\right.$$
If we set $x_0=1$, the first of these recurrences holds for $n=1$ as well as for $n>1$, and we have
$$x_n=1+\sum_{k=0}^{n-1}4^k\;.$$
If this isn’t immediately clear, note that
$$\begin{align*}
x_n&=4^{n-1}+x_{n-1}\\
&=4^{n-1}+4^{n-2}+x_{n-2}\\
&\;\vdots\\
&=4^{n-1}+4^{n-2}+\ldots+4^0+x_0\\
&=1+\sum_{k=0}^{n-1}4^k\;,
\end{align*}$$
and the formula can be proved rigorously by induction on $n$.
Finally, this is just a geometric series, so
$$x_n=1+\frac{4^n-1}{4-1}=1+\frac13(4^n-1)=\frac13(4^n+2)\;.$$
Best Answer
I wrote a little Python program that did the brute force. It found $670$ solutions with $a \le b \le c$ The first was $$a=2019, b= 4074343, c= 16600266807306$$ There were $40$ with $a=2019$ and the last of those was $$a=2019, b=c=8148684$$
As long as nobody laughs too loud, here is the code. I just let $a$ range from $n+1$ to $3n$, compute the range of $b$ to be so that $\frac 1b \le \min(\frac 1a, \frac 1n-\frac 1a)$, compute $c$ using integer division, then see if it comes out even. succ counts the successes. This is Python 2.