From the OP and Zach L's comment to the OP, we can continuously extend $f$ to a function on the extended real numbers by setting $f(-\infty) = 0$ and $f(+\infty) = +\infty$.
Define a sequence $a_n$ of extended real numbers for all natural numbers by $a_0 = -\infty$, and $a_{n+1} = f(a_n)$. Observe that $f$ is a bijective function $[a_n, a_{n+1}] \to [a_{n+1}, a_{n+2}]$.
Let $g$ be the inverse of $f$ (with domain the non-negative extended real numbers)
The sequence $a_n$ is monotonic and increasing, and therefore has a limit $L$ in the extended real numbers. This satisfies
$$\begin{align} L &= \lim_{n \to +\infty} a_n
\\&= \lim_{n \to +\infty} a_{n+1}
\\&= \lim_{n \to +\infty} f(a_n)
\\&= f(\lim_{n \to +\infty} a_n)
\\&= f(L)
\end{align} $$
The OP has already shown that $f$ has no finite fixed points, so therefore $L = +\infty$.
This means the intervals $[a_n, a_{n+1}]$ cover the entire range $[-\infty, +\infty)$.
The fact that $f(f(x)) = x f(x) + 1$ means that the value of $f$ on $[a_{n+1}, a_{n+2}]$ is determined by its values on $[a_n, a_{n+1}]$ (by considering $x \in [a_n, a_{n+1}]$).
Therefore, $f$ is completely determined by its values on $[a_0, a_1] = [-\infty, 0]$.
Conversely, I assert that if you choose any continuous, monotonically increasing function $f_0$ on $[-\infty, 0]$ such that $f_0(-\infty) = 0$ and $0 < f_0(0) < 1$, then
we have a n increasing sequence (converging to $+\infty$) recursively defined by
- $a_0 = -\infty$
- $a_1 = 0$
- $a_2 = f_0(0)$
- $a_{n+2} = a_{n+1} a_n + 1$
and a sequence of invertible functions $f_n : [a_n, a_{n+1}] \to [a_{n+1}, a_{n+2}]$ recursively defined by
- $f_{n+1}(x) = f_n^{-1}(x) x + 1 $
and then the function
$$ f(x) = \begin{cases} f_n(x) & x \in [a_n, a_{n+1}]
\\ +\infty & x = +\infty \end{cases} $$
is (the continuous extension to the extended real numbers of) a solution to the problem.
You can easily check that if we define
$$f_{A} = \frac{ax+b}{cx+d} \quad \text{for} \quad A = \begin{pmatrix}a & b \\ c & d \end{pmatrix}, $$
then $f_{A} \circ f_{B} = f_{AB}$. Thus any matrix $A$ satisfying
$$ A^{2} = k \begin{pmatrix} 1 & 1 \\ 1 & 2 \end{pmatrix}, \quad \text{for some} \ k \neq 0 $$
gives rise to a solution. Mathematica yields two different solutions
$$ f(x) = \frac{1}{x+1} \quad \text{and} \quad f(x) = \frac{2x+1}{x+3}, $$
but I'm not sure if other solutions exist.
Best Answer
HINT: You haven't used the monotonic nondecreasing condition yet.
You have $f(1)=2$ and $f(4)=5$ so $2\le f(2)\le f(3)\le 5$
If you have values up to $f(2n)$ you can use the approach you have for $f(6)$ to find $f(2n+2)$, and then use the monotonic property to find $f(2n+1)$.