Bounds (and range) of a nonlinear difference equation
From experimentation with complex orbits, it looks like the points travel on a curve which is symmetric in $x,y$ and has degree $2$ in $x$. So the curves have equations of the form $a_0 + a_1(x+y) + a_2xy + a_3(x+y)^2 + a_4xy(x+y) + a_5x^2y^2 = 0$.
It is very unlikely that $a_5$ is nonzero because this would give lots of genus $3$ curves (assuming they're not all singular), with a transformation that gives infinitely many rational points on them except when the orbits are periodic, and that can only happen for a countable number of curves. So it would contradict Faltings's theorem. So first I assume $a_5=0$.
Then, since the map $(x,y) \mapsto ((c+y)/x,y)$ preserves them, for a given $y$, the two corresponding $x$ solutions multiply to $c+y$, which means the equation is of the form $b_2(y)x^2 + b_1(y)x + b_0(y)$ where $b_0(y)=(c+y)b_2(y)$, and so $(a_0+a_1y+a_3y^2) = (c+y)(a_3+a_4y) = ca_3+(ca_4+a_3)y+a_4y^2$.
Therefore $a_3=a_4, a_1=(c+1)a_4$ and $a_0=ca_4$, which gives $a_4(c + (c+1)(x+y) + (x+y)^2+xy(x+y)) + a_2xy = 0$,
which means $\frac{c+(c+1)(x+y)+(x+y)^2+xy(x+y)}{xy} = -a_2/a_4$ is constant.
So, the transformation $(x,y) \mapsto (y,(c+y)/x)$ conserves the quantity $\frac{c(1+x+y) + (x+y)(1+x)(1+y))}{xy}$ whenever $xy(c+y) \neq 0$
For almost all fixed values of $c$ and $k$, the curve $\frac{c(1+x+y) + (x+y)(1+x)(1+y))}{xy} = k$ should be an elliptic curve, of which this map is an automorphism. (so most likely a translation of infinite order)
If you start from $(1,1)$, then the curve you obain has the equation $c(1+x+y)+(x+y)(1+x)(1+y) = (3c+8)xy$, or $f(x,y) = c+(c+1)x+(c+1)y+xx-(3c+6)xy+yy+xxy+xyy = 0$.
To get the highest point, you need to intersect it with the curve $\frac{df}{dx} = 0$, so
$c+1+2x-(3c+6)y+2xy+y^2=0$
You get $(2+2y)x = (3c+6)y-(c+1+y^2)$.
Plugging this into the original equation gives you a degree $5$ equation on $y$. $y= -1$ will always be a root because the curve has an horizontal asymptote there, leaving $4$ other real roots, among which you want the second and third largest.
If (like me) you're interested in the dynamics and orbits shapes of this map in the other parts of the plane,
when $k = -1$ the equation factors into $(x+y+1)(c+xy+x+y)$, so a line and an hyperbola, both genus $0$ curves. The map switches one with the other, and doing it twice is doing an automorphism of any of those two curves. Dynamics of automorphisms of genus $0$ curves are well understood. If $c > 3/4$ the two fixpoints are real, so a point gets attracted into a $2$-cycle. If $c < 3/4$ then they are complex, so the map acts more like a "rotation" of the curves. If $c = 3/4$ we get a double (real) fixpoint, so it's more like a "translation".
when $k = -c$ the equation factors into $(x+1)(y+1)(c+x+y)$, so three lines. The line $y=-1$ is sent onto the line $x=-1$ then to $x+y=-c$ then to $y=-1$ again. In this case, the discriminant is $(c-2)²-3$, so the behaviours change at $c= 2\pm \sqrt 3$
when $k = \infty$ the equation becomes $xy=0$ (or $xyz = 0$ if you're working in the projective plane), so three lines. The line $y=0$ is sent to the line $x=0$ then to the blowup of a point at infinity, then to the line at infinity $(z=0)$ then to another blowup then back to $y=0$. This long trip sends $(x,0)$ to $(x/c,0)$. Thus we get an attractive $5$-cycle most of the time.
There are another two values of $k$ (the roots of $8-k-13c+10ck-ck^2+16c^2$) where the curve is singular (there is a node where the small loops disappear), and so the smooth locus of the curve has genus $0$ that is fixed by the transformation.
For any other values of $k$, the curve is an elliptic curve (with one or two real components), and the map is some infinite order translation in the generic case.
If you compute the $j$ invariants, those special values of $k$ correspond to poles, and their order is the number of genus $0$ curves the transformation cycles through. The map $k \mapsto j(C_k)$ is then a rational function of degree $5+3+2+1+1 = 12$.
All the curves go through the points $(-1,0),(-1,c),(0,-1),(0,-c)$, have common vertical and horizontal asymptotes (so they all go to the same point on the two blowups I mentioned before), plus they each have their own oblique asymptote, so that's another common point on the line at infinity.
If you blow up the projective plane at those $4$ finite points, those $3$ points at infinity, then blowup again at the two common asymptotes, the original map should maybe become a nice automorphism of a smooth elliptic surface
As for periodicity, I expect the set of $c$ where this happens is dense (and those values of $c$ are all algebraic). One would need to investigate how the $j$-invariant and the translation varies while $k$ and $c$ varies.
Computational evidence suggest the following conjectures.
For $1\ne c>0$ a sequence $\{z_n\}=\{(x_n,x_{n+1})\}$ is a dense set of a smooth curve which bounds a convex shape. When $c$ tends to infinity this shape tends to an isosceles right angled triangle with sides of order $c$, when $c$ tends to zero the shape tends to an ellipse, see the graphs.
$c=2^{-2}, 2^{-1},\dots, 2^7$ (I remark that the diameter of the shape is not a monotone function of $c$):
$c=2^{8}, \dots, 2^{11}$:
$c=2^{-3}, \dots, 2^{-22}$:
Everybody feel free to find the equations of the curves and to prove the conjectures, I'll try to do this too. In particular, if these equations are polynomial then their coefficients can be found by solving the respective system of linear equations.
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)$.