If $a_{1}=1$ and $a_{n+1}=1+\frac{n}{a_{n}}$, then $a_{n}=\sqrt{n}+\frac{1}{2}-\frac{1}{8\sqrt{n}}+o\left(\frac{1}{\sqrt{n}}\right)$ let sequence $\{a_{n}\}$ such $$a_{1}=1,a_{n+1}=1+\dfrac{n}{a_{n}}$$
show that:
$$a_{n}=\sqrt{n}+\dfrac{1}{2}-\dfrac{1}{8\sqrt{n}}+o\left(\dfrac{1}{\sqrt{n}}\right)?$$
 A: It will be simpler in what follows to use two steps at once, so we can start by computing the relationship $a_{n+2} = \frac{a_n(n+2) + n}{a_n + n}$, and letting $b_n = a_n / \sqrt n$ we get $b_{n+2} = \frac {b_n (1 + 2/n) + 1/\sqrt n}{b_n\sqrt{1/n+2/n^2} + \sqrt{1+2/n}}$.
Let's denote $\begin{bmatrix} a & b \\ c & d \end{bmatrix} x = \frac {ax+b}{cx+d}$.
If we have a relation $b_{n+2} = \begin{bmatrix} 1+A & B \\ C & 1+D \end{bmatrix} b_n$ where each $A,B,C,D$ is an $O(1/\sqrt n)$ and we look then at some $c_n = (b_n-k)\sqrt n$, we obtain on the sequence $c_n$ the same kind of relation (so with dominant term still begin $I_2$) with
$$A' = (1+A-kC)\sqrt{1+2/n}-1, \\B' = (B-kD+kA-k^2C)\sqrt{n+2}, \\
C' = C/\sqrt n, \\ D' = D+kC$$.
Right now, $C$ is still as big as everyone else, so to make $B'$ not an order of magnitude bigger than the rest, we have to solve a degree $2$ equation in $k$, which gives $k= \pm 1$.
If we define $c_n = (b_n-1)\sqrt n$, we get $A' = -n^{-1/2} + O(n^{-1}), D' = +n^{-1/2} + O(n^{-1}), C' = O(n^{-1}), B' = O(n^{-1/2})$.
From now $C'$ will continue to plummet, the dominant terms of $A'$ and $D'$ will not be able to change, and so to make the next $B$ into an $O(n^{-1/2})$ term, we will have to pick $k = $ the constant term of $B/(A-D)$.
What this tells us is that there are exactly two asymptotic developments for $a_n$ whose error terms at each level obey a recurrent relation of the form $r_{n+2} = [I_2 + O(n^{-1/2})] r_n$.
Since all of this can be done algebraically, those two power series are $\sqrt n$ times the two solutions in $\Bbb Q[[n^{-1/2}]]$ to the original equation on $a_n/\sqrt n$ (and you go from one to the other by switching the sign of $\sqrt n$).
Now we can look at $u_n = ((((a_n/\sqrt n - 1)\sqrt n - \frac 12)\sqrt n + \frac 18)\sqrt n + \frac 18)\sqrt n$.
It satisfies $u_{n+2} = \begin{bmatrix}1 - n^{-1/2} + O(n^{-1}) && \frac 7 {64}n^{-1} + O(n^{-3/2}) \\ O(n^{-5/2}) && 1 + n^{-1/2} + O(n^{-1}) \end{bmatrix} u_n$,
(and each constant in those $O$ can be effectively computed)  
This means that if you have $u_n$ close enough to $0$, then $u_{n+2} = (1-2n^{-1/2})u_n + O(n^{-1})$. So it should mean that for each bound $\delta$ there is an $n_0$ such that $\forall n > n_0, |u_n < \delta | \implies |u_{n+2} < \delta|$.
If you effectively compute $n_0$ for, say, $\delta = 1$, this allows you to provably check that a sequence $u_n$ stays bounded by finding a term $u_n$ with $n$ large enough that is below the $\delta$. If the asymptotic development is valid, a finite computation will be able to prove it.

I would love to see a proof that the asymptotic development is valid for every starting value of $a_n$ (except the one value that will stay away infinitely from the asymptotic formula and will instead obey the other asymptotic development)
A: I will present two approaches to this question. The first is less advanced and easier but also yields a less optimal result.

First Method
Consider
 the auxiliary sequence $b_n =\sqrt{n+\frac{1}{4}}+\frac{1}{2}$. This is
 the unique positive solution to the $P_n(X)=0$ where $P_n(X)=X^2-X-n$.

${\bf Lemma~1.}$ For every $n\geq1$, we have
  $b_{n+1}\geq 1+\dfrac{n}{b_{n-1}}$.

Proof. Indeed, let $u=1+\frac{n}{b_{n-1}}$ then
$$u^2-u=\frac{(n+b_{n-1})n}{b_{n-1}^2}=\frac{(n+b_{n-1})n}{b_{n-1}+n-1}=
n+\frac{n}{b_{n-1}+n-1}\leq n+1$$
since $b_{n-1}\geq 1$. This shows that $P_{n+1}(u)\leq0$ and 
consequently $u\leq b_{n+1}$ as desired.$\qquad \square$

${\bf Lemma~2.}$ For every $n\geq1$, we have
  $b_{n-1}\leq a_n \leq b_n$.

Proof. This is now an easy induction. Clearly true for $n=1$, and if 
it is true for some $n$ then
$$
b_n=1+\frac{n}{b_n}\leq 1+\frac{n}{a_n}\leq 1+\frac{n}{b_{n-1}}\leq 
b_{n+1}.
$$
and the result follows.$\qquad \square$
It follows that
$$
\sqrt{n-\frac{3}{4}}-\sqrt{n}\leq a_n-\sqrt{n}-\frac{1}{2}\leq \sqrt{n+\frac{1}{4}}-\sqrt{n}
$$
Thus
$$\sqrt{n}\left\vert a_n-\sqrt{n}-\frac{1}{2}\right\vert\leq \frac{3}{4(1+\sqrt{1-3/(4n)})}\leq \frac{1}{2}$$
So, we have proved that, for every $n\geq 1$ we have
$$
a_n= \sqrt{n}+\frac{1}{2}+{\cal O} \left(\frac{1}{\sqrt{n}}\right)
$$
This is not the desired expansion but it has the merit to be easy to prove.

Second Method
The general reference in this part is the book of D. Knuth. "The art of Computer programming, Vol III, second edition, pp.63--65".
Let $I(n)$ be the number of involutions  in the symmetric group $S_n$, ($i.e.$ $\sigma\in S_n$ such that $\sigma^2=I$). It is well-known that $I(n)$ can be calculated inductively by
$$
I(0)=I(1)=1,\qquad I(n+1)=I(n)+nI(n-1) 
$$
This shows that our sequence $\{a_n\}$ is related to these numbers by the formula
$$a_n=\frac{I(n+1)}{I(n)}.$$
So, we can use what we know about these numbers, In particular, the following asymptotic expansion, from Knuth's book:
$$
I(n+1)=\frac{1}{\sqrt{2}} n^{n/2}e^{-n/2+\sqrt{n}-1/4}\left(1+\frac{7}{24\sqrt{n}}+{\cal O}\left(\frac{1}{n^{3/4}}\right)\right).
$$
Now, it is a "simple" matter to conclude from this that $a_n$ has the desired asymptotic expansion.
Edit: In fact the term $\mathcal{O}(n^{-3/4})$ effectively destroys the asymptotic expansion as mercio noted, But in fact we have 
$$
I(n+1)=\frac{1}{\sqrt{2}} n^{n/2}e^{-n/2+\sqrt{n}-1/4}\left(1+\frac{7}{24\sqrt{n}}+{\cal O}\left(\frac{1}{n}\right)\right).
$$
This is shown by WIMP AND ZEILBERGER in their paper ``Resurrecting the Asymptotics
of Linear Recurrences'', that can be found here.
