How to prove this inequality $x^2_{n}\le\frac{8}{3}$

Very interesting problem. A possible approach to the solution is through the theory of dynamical systems considering it as a perturbation of a conservative system. I guess that the solution below might be improved to provide a slightly better constant.

Setting $y_n=x_{n-1}$ the given equation is equivalent to $$ \left\{ \begin{array}{l} y_{n+1} = x_n\\ x_{n+1} = \left(1+\frac{1}{n}\right)x_n-y_n \end{array} \right. $$

Consider the following energy function for the system $$ E_n = (x_n + y_n)^2+3(x_n-y_n)^2. $$ Using the definition of $x_{n+1}$ and $y_{n+1}$ above, a somehow lengthy, but straightforward computation gives $$ E_{n+1} = E_n + \frac{4}{n}\left(1+\frac{1}{n}\right)x_n^2 - \frac{8}{n}x_ny_n = E_n + \frac{4 x_n^2}{n^2} + \frac{4}{n} \left(x_n (x_n-y_n) -x_ny_n\right). $$ Recalling that $x_n-y_n = x_{n+1}-\frac{1}{n}x_n$, this implies that $$ E_{n+1} = E_n + \frac{4}{n} (x_{n+1}x_n - x_n y_n) $$ and hence $$ \begin{align} E_{n+1} &= E_2 + \sum_{i=2}^n \frac{4}{i} (x_{i+1} y_{i+1} -x_iy_i) = E_2 + \frac{4}{n} x_{n+1} y_{n+1}+\sum_{i=3}^{n} \left( \frac{4}{i-1} - \frac{4}{i} \right)x_{i}y_{i} - 2x_2y_2\\ &=4+\frac{4}{n} x_{n+1} y_{n+1} + \sum_{i=3}^{n} \frac{4}{i(i-1)} x_{i}y_{i} \end{align} $$ where in the last passage we used $x_2 = 1$ and $y_2 = 0$.

Observing that $E_{n} \geq (x_n+y_n)^2+(x_n-y_n)^2 = 2(x_n^2+y_n^2)$ we deduce $$ 2(x_{n+1}^2 + y_{n+1}^2) \leq E_{n+1} \leq 4 + \frac{2}{n}(x_{n+1}^2 + y_{n+1}^2)+\sum_{i=3}^{n} \frac{2(x_i^2 + y_i^2)}{i(i-1)} $$ or $$ x_{n+1}^2 + y_{n+1}^2 \leq 4 +\sum_{i=3}^{n} \frac{2(x_i^2 + y_i^2)}{i(i-1)} $$

We now want to apply Gronwall's discrete lemma to this equation (see here). Recalling $$ \sum_{i=2}^{\infty}\frac{2}{i(i+1)} = 1 $$ we prove $x_n^2 \leq 4e$ which shows that the desired sequence is bounded despite with a worse constant than asked.

P.S.: Numerical inspection in the $(x,y)$ plane above suggests the optimal bound to be $x_n^2 \leq \frac{9}{4}$

P.P.S.: I was asked to add some insights on how the argument was developed. The techniques used are quite common in the analysis of ODEs. A good, undergraduate, introduction to some of the ideas below can be found in "Nonlinear Dynamics and Chaos" by Strogatz (link). I am not aware of thorough accessible works specifically dedicated to multidimensional discrete dynamical systems (which often arise though as Poincarè maps of higher dimensional systems of ODEs) since already the one-dimensional case gives rise to several tough problems.

Concerning the main steps of the argument:

• The transformation into a first order system is a common trick.
• Direct inspection of the "unperturbed" system (i.e. the one without the $1/n$ term) shows that all sequences generated by arbitrary $x_1$ and $x_2$ are 6-periodic. Hence the idea to look for an almost conserved quantity for the given case. Some educated guessing supported by simulation led to the definition of the "energy" $E_n$ above which is conserved up to $O(\frac{1}{n})$ terms.
• Gronwall's technique is finally standard tool in the field. However, in this case some pre-work was needed on the energy estimate to leave only converging terms (hence the somehow painful rearrangement and index shift in the latest steps).

Just so this question is fully answered (i.e. getting the bound $8/3$) here I add a condensed proof of $x_n^2 \leq \frac{8}{3}$ from the answer given here as Community Wiki.

---

Define $$E_n = (x_n - x_{n-1})^2 + F_n\left(x_n + x_{n-1}\right)^2$$ with $F_n = \frac{(n-1)}{3n+1}$. The recurrence relation gives us $$E_{n+1} - E_n = \frac{(F_{n+1}-F_n)}{F_n}\cdot F_n\left(x_n + x_{n+1}\right)^2 \geq 0$$ since $F_n$ is a positive and increasing function. This shows that $E_n$ is increasing. We also see that $$E_{n+1} - E_n \leq \frac{(F_{n+1}-F_n)}{F_n}\cdot E_n \implies \frac{E_{n+1}}{F_{n+1}} \leq \frac{E_n}{F_n}$$ so $\frac{E_n}{F_n}$ is decreasing and since $F_n$ converges it follows that $E_n$ converges. To bound $x_n$ from $E_n$ note that the definition of $E_n$ (the equation for an ellipse) allows us to write $x_n - x_{n-1} = \sqrt{E_n}\cos(\phi_n)$ and $x_n + x_{n-1} = \sqrt{\frac{E_n}{F_n}}\sin(\phi_n)$ for some angle $\phi_n$. Thus

$$x_n^2 = \frac{E_n}{4}[\cos(\phi_n) + \frac{1}{\sqrt{F_n}}\sin(\phi_n)]^2$$ Maximizing the right hand side over $\phi_n$ gives us $$x_n^2 \leq \frac{E_n}{4}\left[1 + \frac{1}{F_n}\right] = E_n\cdot \frac{n}{n-1}$$ Since $\frac{E_n}{F_n}$ is decreasing we have $E_n \leq F_n \frac{E_2}{F_2} = 8F_n$ and $$x_n^2 \leq \frac{8n}{3n+1} \leq \frac{8}{3}$$ It is possible to improve the constant with a bit of numerical computation (using $n=5$ instead of $n=2$ when bounding $E_n$) to the best possible one $\frac{9}{4}$ (for which we have equality for $n=3$). The best possible asymptotic bound would be slightly smaller than $2$ (best possible in the sense that $x_n^2 \approx 2$ would occur infinitely many times).

(Too long for a comment)

Let $S_0 = 0, \ S_n = \sum_{k=1}^n x_k.$ Then applying summation by parts to the formula

$$ x_{n+1}=\sum_{k=1}^{n}\dfrac{x_{k}}{k}-\sum_{k=1}^{n-1}x_{k}$$

yields the equation (for $n\geq 2$)

$$ x_{n+1} = 1 + \frac{S_n}{n} - S_{n-1} + \sum_{k=1}^{n-1} \frac{S_k}{k(k+1)}$$

Adding $S_n$ to both sides of this also gives

$$S_{n+1} = 1 + \frac{S_n}{n} + x_n + \sum_{k=1}^{n-1} \frac{S_k}{k(k+1)}$$

I attempted to use these equations together with an inductive hypothesis of the form $S_{n-1} \in [a,b] , x_n \in [c,d]$ to show that $x_{n+1}, S_{n+1}$ must also lie in the same bounds. While you can get extremely close, it never quite works out and I've convinced myself that this method can not succeed no matter with choices we make for $a,b,c,d.$

However, using an established a bound $|x_n| \leq M$ (such as the one Stefano proved above) then we can prove inductively that $S_n$ is also bounded.

Conjecture - The sequence $S_n$ approaches a limit $L,$ and $$L = 1 + \sum_{n=1}^{\infty} \frac{S_n}{n(n+1)} \approx 1.953053682$$ A corollary of this is that $x_n \to 0,$ which it appears to do so very very slowly.