<s>[INST] Bounds (and range) of a nonlinear difference equation I'm interested in the following set of nonlinear difference equations:
$$x_{n+1} = \frac{c + x_n}{x_{n-1}},\; x_1 = x_0 = 1 \qquad \textrm{for } c > 0$$ 
For $c=1$ the sequence is periodic with period 5 and range {1,2,3}. For other values of $c$ the sequence appears to be oscillatory. For small $c$, the supremum seems to be slightly more than $1+2c$, while for larger $c$ it seems to approach $3c$.
Is there a simple formula for the bounds of the sequence in this general case? What about the range? Are the sequences periodic? Are they dense within the bounds? 

 [/INST] This is an partial answer. Please consider this as a supplement to merico's excellent answer.
Introduce auxillary sequences $(y_n), (z_n), (w_n)$ by $y_n = x_{n+1}, z_n = x_{n+2}, w_n = x_{n+3}$.
Rewrite the defining recurrence relation of $(x_n)$ to a more symmetric form:
$$x_{n+1} = \frac{c+x_n}{x_{n-1}}\quad \iff\quad x_{n+1}x_{n-1} = c + x_n$$
In terms of the auxillary sequences, this becomes
$$x_n z_n = y_n + c\quad\text{ and }\quad y_n w_n = z_n + c.$$ 
As a result,
$$
(1+x_n)z_n = z_n + z_n x_n = z_n + (y_n+c) = y_n + (z_n + c) = y_n + y_nw_n = y_n(1+w_n)
$$
Multiply both side by $\frac{(1+y_n)(1+z_n)}{y_nz_n}$, we find
$$\frac{(1+x_n)(1+y_n)(1+z_n)}{y_n}
= \frac{(1+y_n)(1+z_n)(1+w_n)}{z_n}
= \frac{(1+x_{n+1})(1+y_{n+1})(1+z_{n+1})}{y_{n+1}}$$
This means following expression
$$\frac{(1+x_n)(1+y_n)(1+z_n)}{y_n} = \frac{(1+x_n)(1+y_n)(x_n+y_n+c)}{x_ny_n}$$
is independent of $n$ and they are invariant of the system.
Up to an additive offset $c$, the second form is the invariant found by merico.
Since $x_1 = y_1 = 1, z_1 = c+1$, the sequence of points $(x_n, y_n, z_n)$ lies on intersection of following two surfaces in $\mathbb{R}^3$:
$$(1+x)(1+y)(1+z) = 4(c+2)y\quad\text{ and } xz = y + c$$
Parametrize the curves within above interesection by some parameter $s$ and taking logarithm derivatives, we obtain
$$\frac{x'}{1+x} + \frac{y'}{1+y} + \frac{z'}{1+z} = \frac{y'}{y}\quad\text{ and }\quad
\frac{x'}{1+x} + \frac{z'}{1+z} = \frac{y'}{y+c}
$$
At those $s$ where where $y$ is extremal, $y' = 0$. 
Above equations implies $x = z$ and hence $x^2 = y + c$. Furthermore,
$$(1+x)^2(1+y) = 4(c+2)y
\iff (1+2x+y+c)(1+y) = 4(c+2)y\\
{\large\Downarrow}\\
x = -\frac{y^2 - 3(c+2)y + (c+1)}{2(y+1)}
$$
Substitute this back into $x^2 = y + c$, we find the extremal values of $y$
are roots of a quartic polynomial.
$$(y^2 - 3(c+2)y + (c+1))^2 - 4(y+1)^2(y+c) = 0 \tag{*1}$$
For $c > 0$, this polynomial has $4$ distinct real roots: $0 \le \lambda_1(c) < \lambda_2(c) < \lambda_3(c) < \lambda_4(c)$. 
Start from $x_1 = x_2 = 1$, it is easy to see all $x_n, y_n$ generated are positive. This means the extremal values of $y$ also need to satisfy:
$$y^2 - 3(c+2) + (c+1) \le 0
\iff \frac32(c+2) - \Delta(c) \le y \le \frac32(c+2) + \Delta(c)\tag{*2}
$$
where $\Delta(c) = \frac12\sqrt{9c^2+32(c+1)}$.
If one make an implicit plot of the roots of $(*1)$. One find in general,
$$\lambda_1(c) < \frac32(c+2) - \Delta(c) < \lambda_2(c) < \lambda_3(c) < \frac32(c+2) + \Delta(c) < \lambda_4(c)$$
From this, we can deduce in general, $x_{n+1} = y_n$ are bounded between
the two middle two roots $\lambda_2(c), \lambda_3(c)$ of equation $(*1)$.
</s>