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? 


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)$.
