Asymptotic analysis of $x_{n+1} = \frac{x_n}{n^2} + \frac{n^2}{x_n} + 2$ 
Problem: Let $x_1 = 1$ and $x_{n+1} = \frac{x_n}{n^2} + \frac{n^2}{x_n} + 2, \ n\ge 1$.
Find the third term in the asymptotic expansion of $x_n$.

I have posted it in MSE six months ago without solution for the third term
https://math.stackexchange.com/questions/3801405/the-limit-and-asymptotic-analysis-of-a-n2-n-from-a-n1-fraca-nn.
We have $\lim_{n\to \infty} (x_n - n) = \frac{1}{2}$ (see [1]; I also give a solution with the help of computer in the link above).
So the first two terms in the asymptotic expansion of $x_n$ are $x_n \sim n + \frac{1}{2}$.
Edit: In [1], the authors proved that $\frac{1}{4n-2} \le x_n - n - \frac{1}{2} \le \frac{2}{2n-3}$ for all $n\ge 3$.
For the third term, @Diger in MSE said $x_n \sim n + \frac{1}{2} + \frac{5}{8n}$ (see @Diger's answer in the link above).
However, I did some numerical experiment which does not support this result.
I am not convinced of the numerical evidence due to finite precision arithmetic.
I hope to prove or disprove it analytically.
Numerical Experiment: If $x_n \sim n + \frac{1}{2} + \frac{5}{8n}$,
then it should hold $16n(x_{2n} - 2n - \frac{1}{2}) \approx 5$ and $8(2n+1)(x_{2n+1} - (2n+1) - \frac{1}{2}) \approx 5$
for large $n$. When $n=1500$, I use Maple to get
$16n(x_{2n} - 2n - \frac{1}{2}) \approx 4.368$ and $8(2n+1)(x_{2n+1} - (2n+1) - \frac{1}{2}) \approx 5.642$.
When $n$ is larger (e.g., $n=10000$), the numerical result seems unreliable.
I ${\color{blue}{\textbf{GUESS}}}$ that
$$x_{2n} \sim 2n + \frac{1}{2} + \frac{q_1}{2n},$$
$$x_{2n+1} \sim (2n+1) + \frac{1}{2} + \frac{q_2}{2n+1}$$
where $q_1 + q_2 = \frac{5}{4}$ and $q_1 \ne q_2$ (if $q_1 = q_2$, then it is $x_n \sim n + \frac{1}{2} + \frac{5}{8n}$).
(Some numerical experiment shows $q_1 \approx \frac{61}{112}, q_2 \approx \frac{79}{112}$. But I am not convinced of it.)
Edit: I give more analysis for my guess as an answer.
Any comments and solutions are welcome and appreciated.
Reference
[1] Yuming Chen, Olaf Krafft and Martin Schaefer, “Variation of a Ukrainian Olympiad Problem: 10982”,
The American Mathematical Monthly, Vol. 111, No. 7 (Aug. - Sep., 2004), pp. 631-632
 A: Consider the substitutions
\begin{equation*}
    x_n=n+1/2+y_n/n,\quad y_n=u_n+5/8. 
\end{equation*}
Then $u_1=-9/8$ and
\begin{equation*}
    u_{n+1}=f_n(u_n)
\end{equation*}
for $n\ge1$, where
\begin{equation*}
    f_n(u):=\frac{-64 n^4 u-8 n^3 (4 u-13)+n^2 (56 u+115)+n (96 u+76)+4 (8 u+5)}{8 n^2 \left(8 n^2+4 n+8 u+5\right)}. 
\end{equation*}
Define $c_n(u)$ by the identity
\begin{equation*}
    f_n(u)=-u+\frac{13}{8n}+\frac{c_n(u)}{n^2}, 
\end{equation*}
so that
\begin{equation*}
    c_n(u)=\frac{n^2 \left(64 u^2+96 u+63\right)+n (11-8 u)+4 (8 u+5)}{8 \left(8 n^2+4 n+8 u+5\right)}. 
\end{equation*}
Then for $n\ge1$
\begin{equation*}
    u_{n+1}+u_n=\frac{13}{8n}+\frac{c_n(u)}{n^2} \tag{1}
\end{equation*}
and for $n\ge2$
\begin{equation*}
    u_{n+1}=f_n(f_{n-1}(u_{n-1}))=u_{n-1}-\frac{13}{8n(n-1)}+\frac{c_n(u_n)}{n^2}
    -\frac{c_{n-1}(u_{n-1})}{(n-1)^2}. \tag{2}
\end{equation*}
Note that
\begin{equation*}
    u_{101}=-0.54\ldots,\quad u_{102}=0.56\ldots,  \tag{3}
\end{equation*}
and
\begin{equation*}
    0\le c_n(u)\le3  
\end{equation*}
if $n\ge10$ and $u\in[-6/10,8/10]$. Therefore and because for natural $m\ge102$ we have
\begin{equation*}
    \sum_{n=m}^\infty\Big(\frac{13}{8n(n-1)}+\frac3{(n-1)^2}\Big)<\frac5{m-2}\le0.05,
\end{equation*}
it follows from (2) and (3) by induction that for all $n\ge101$ we have $u_n\in[-6/10,8/10]$ and hence $0\le c_n(u_n)\le3$. So, again by (2), the sequences $(u_{2m})$ and $(u_{2m+1})$ are Cauchy-convergent and hence convergent. Moreover, by (1), $u_{n+1}+u_n\to0$.
Thus, indeed
\begin{equation*}
    y_{n+1}+y_n\to5/4, 
\end{equation*}
and the sequences $(y_{2m})$ and $(y_{2m+1})$ are convergent. (The limits of these two sequences can in principle be found numerically with any degree of accuracy -- controlled by (2), say.)
