Empy2 pointed out the error in your proof. You cannot deduce that $f\big(f(x)f(0)\big)+f(x)=f(0)$ for all $x$ implies $f\big(xf(0)\big)+x=f(0)$ for all $x$. This only works only for $x$ in the range of $f$. I am presenting a different solution.
If $f(0)=0$, as you showed, we have
$$f(0)+f(x)=f\big(f(x)f(0)\big)+f(x+0)=f(x\cdot 0)=f(0)$$
so $f(x)=0$ for all $x$. We assume that $f(0)\ne 0$ from now on. We claim that $f(c)=0$ implies that $c=1$.
Suppose that $f(c)=0$ for some $c\ne 1$. Then, we have
$$f\Biggl(f(c)f\left(\frac{c}{c-1}\right)\Biggr)+f\left(c+\frac{c}{c-1}\right)=f\Biggl(c\left(\frac{c}{c-1}\right)\Biggr).$$
That is,
$$f(0)+f\left(\frac{c^2}{c-1}\right)=f\left(\frac{c^2}{c-1}\right),$$
or $f(0)=0$, contradicting the assumption that $f(0)\ne 0$. Therefore, the hypothesis that $c\ne1$ is false.
Now, we note that for every $x\ne1$,
$$f\Biggl(f(x)f\left(\frac{x}{x-1}\right)\Biggr)+f\left(x+\frac{x}{x-1}\right)=f\Biggl(x\left(\frac{x}{x-1}\right)\Biggr),$$
so
$$f\Biggl(f(x)f\left(\frac{x}{x-1}\right)\Biggr)+f\left(\frac{x^2}{x-1}\right)=f\left(\frac{x^2}{x-1}\right).$$
So,
$$f\Biggl(f(x)f\left(\frac{x}{x-1}\right)\Biggr)=0$$
and we then conclude from our claim above that
$$f\left(\frac{x}{x-1}\right)=\frac{1}{f(x)}\tag{1}$$
for every $x\ne 1$. In particular, this shows that $f(1)=0$, since there does exist $c$ such that $f(c)=0$.
From (1), we get $f(0)=\pm 1$. Observe that $f$ is a solution if and only if $-f$ is also a solution. We can without loss of generality assume that $f(0)=-1$. Plugging $y\mapsto 1$ in the original functional equation, we have
$$-1+f(x+1)=f(0)+f(x+1)=f\big(f(x)f(1)\big)+f(x+1)=f(x\cdot1)=f(x),$$
or
$$f(x+1)=f(x)+1.\tag{2}$$
By induction, $f(x+n)=f(x)+n$ for all integers $n$.
Here, we claim that $f$ is injective. Suppose that $f(u)=f(v)$ for some $u,v\in\mathbb{R}$. Pick a positive integer $n$ so large that $4(u+n)< (v+n+1)^2$. Hence, there are two distinct $a,b\in\mathbb{R}$ such that $u+n=ab$ and $v+n+1=a+b$ (since $a$ and $b$ are the roots of the polynomial $t^2-(v+n+1)t+(u+n)$, which have two distinct real roots). Therefore,
$$f\big(f(a)f(b)\big)+f(a+b)=f(ab)=f(u+n)=f(u)+n.$$
But $f(a+b)=f(v+n+1)=f(v)+n+1=f(u)+n+1$. That is, by (2), we have $$f\big(f(a)f(b)+1\big)=f\big(f(a)f(b)\big)+1=0.$$ This means $f(a)f(b)+1=1$, or $f(a)f(b)=0$. Consequently, $f(a)=0$ or $f(b)=0$, which means $a=1$ or $b=1$. Without loss of generality, $b=1$, so $u+n=ab=a$ and $v+n+1=a+b=a+1$. That is, $u=a-n=v$.
Now, substitute $y\mapsto 1-x$ in the original functional equation. We then have
$$f\big(f(x)f(1-x)\big)=f\big(f(x)f(1-x)\big)+f(1)=f\big(x(1-x)\big).$$
Thus, by injectivity,
$$f(x)f(1-x)=x(1-x).$$
Then, we take $y\mapsto -x$ in the original functional equation to get
$$f\big(f(x)f(-x)\big)-1=f\big(f(x)f(-x)\big)+f(0)=f\big(x(-x)\big).$$
Thus,
$$f\big(f(x)f(-x)\big)=f(-x^2)+1=f(-x^2+1)$$
by (2). By injectivity,
$$f(x)f(-x)=-x^2+1.$$
From (2), we also have
\begin{align}f(x)&=f(x)\Big(\big(f(-x)+1\big)-f(-x)\Big)\\&=f(x)\big(f(1-x)-f(-x)\big)=f(x)f(1-x)-f(x)f(-x),\end{align}
so
$$f(x)=x(1-x)-(-x^2+1)=x-1.$$
It is easy to see that $f(x)=x-1$ is indeed a solution. Therefore, all solutions are $f(x)=0$, $f(x)=x-1$, and $f(x)=1-x$.
You can show that the functions $ f : \mathbb R \to \mathbb R $ satisfying
$$ f \big( x + f ( y ) \big) = y + f ( x + 1 ) \tag 0 \label 0 $$
for all $ x , y \in \mathbb R $ are exactly those of the form $ f ( x ) = A ( x ) + 1 $, where $ A $ is an additive involution, i.e. $ A ( x + y ) = A ( x ) + A ( y ) $ and $ A \big( A ( x ) \big) = x $ for all $ x , y \in \mathbb R $. If further regularity conditions like continuity, local boundedness, integrability or measurability are assumed for $ f $, then $ A $ will be regular, too. The only regular additive functions are those of the form $ A ( x ) = c x $ for some constant $ c \in \mathbb R $ (see here). Thus the only regular additive involutions are $ A ( x ) = x $ and $ A ( x ) = - x $, as we must have $ A \big( A ( 1 ) \big) = c ^ 2 = 1 $. Therefore, the only regular solutions to \eqref{0} are $ f ( x ) = x + 1 $ and $ f ( x ) = - x + 1 $. Without any further conditions on $ f $, one can use the axiom of choice to show that there are nonregular solutions, too (see Example 5.6 in this PDF file).
It's straightforward to check that if $ f $ is of the form $ f ( x ) = A ( x ) + 1 $ for some additive involution $ A $, then it satisfies \eqref{0}. We try to prove the converse.
Letting $ x = 0 $ in \eqref{0} we get
$$ f \big( f ( y ) \big) = y + f ( 1 ) \text . \tag 1 \label 1 $$
In particular, \eqref{1} shows that if $ f ( x ) = f ( y ) $ then $ x = y $. Letting $ y = 0 $ in \eqref{1} shows that $ f \big( f ( 0 ) \big) = f ( 1 ) $, and thus by injectivity, $ f ( 0 ) = 1 $. Hence, putting $ x = - 1 $ in \eqref{0} we have
$$ f \big( f ( y ) - 1 \big) = y + 1 \text . \tag 2 \label 2 $$
\eqref{2} shows that if we substitute $ x - 1 $ for $ x $ and $ f ( y ) - 1 $ for $ y $ in \eqref{0}, we get
$$ f ( x + y ) = f ( x ) + f ( y ) - 1 \text . \tag 3 \label 3 $$
Now, if we define $ A : \mathbb R \to \mathbb R $ with $ A ( x ) = f ( x ) - 1 $, then by \eqref{3} we can see that $ A $ is additive, and by \eqref{2} we can see that $ A $ is an involution. This proves what was desired.
Best Answer
There is indeed no such function
Suppose such $f : (1, \infty)\rightarrow(1, \infty)$ exists. Observe firstly that the dimensionality of the constraint on $f$ causes $f$ to be determined by its value at any one point in the following way. Let's consider, arbitrarily, $f(2)$. Let $f(2) = a$, then we have by the constraint for $y = 2$:
$$x^{(a^x)} = 2^{(f(x)^2)}$$ $$a^x \log_2 x = f(x)^2$$ $$f(x) = \sqrt{a^x \log_2 x}$$
Now we see if this holds up to scrutiny elsewhere along the constraint. Note the constraint is trivial if $x = y$, and of course it will hold for $y = 2$ (and by symmetry, $x = 2$) by design, so let's try, say $x = 4$ and $y = 8$. The constraint is then: $$4^{\left(\sqrt{a^8 \log_2 8}\right)^4} = 8^{\left(\sqrt{a^4 \log_2 4}\right)^8}$$ $$4^{\left(\sqrt{3a^8}\right)^4} = 8^{\left(\sqrt{2a^4}\right)^8}$$ $$4^{(9a^{16})} = 8^{(16a^{16})}$$ $$9a^{16}\log_2 4 = 16a^{16}\log_2 8$$ $$18a^{16} = 48a^{16}$$
Clearly this has the sole solution $a = 0$. Yet, by supposition, $a$ is in the range of $f$, but this contradicts $f$ being into $(1, \infty)$, thus no such $f$ exists.