Fdtxjr

A closed form I doubt there is. But asymptotics are easy:
$$
a_n+1^2=a_n^2+2+1/a_n^2,
$$
hence, for every $nge1$,
$$
a_n^2=2n+1+sum_k=0^n-1frac1a_k^2.qquadqquadqquadqquad (*)
$$
This shows that $a_n^2ge2n+2$ for every $nge1$, for example $a_100gesqrt202>10sqrt2>14$. In particular, $a_nto+infty$. Plugging this into $(*)$ yields $a_n^2=2n+1+o(n)$ hence
$$
sqrt2nle a_nlesqrt2n+o(sqrtn).
$$
At this point, we know that $a_n^2ge2n+2$ for every $nge1$. Using $(*)$ again, one sees that, for every $nge1$,
$$
a_n^2le2n+2+sum_k=1^n-1frac12+2kle2n+2+frac12log(n).
$$
Which already shows that $$14.2<a_100<14.3$$
Plugging this upper bound of $a_n^2$ into $(*)$ would yield a refined lower bound of $a_n^2$. And one could then plug this refined lower bound into $(*)$ again to get a refined upper bound. And so on, back and forth between better and better upper bounds and better and better lower bounds. (No more asymptotics here.)

To elaborate on how I got my answer: I started with @Did's $a(n) approx sqrt2n$ and looked for a next term. $a(n) = sqrt2n$ would make $ a(n+1) - (a(n) + a(n)^-1) = sqrt 2,n+2-sqrt 2sqrt n-1/2,frac sqrt 2sqrt n = - fracsqrt28 n^-3/2 + O(n^-5/2)$. With $a(n) = sqrt2n + c n^-1/2$ I don't get a change in the $n^-3/2$ term, so I tried $a(n) = sqrt2n + c ln(n) n^-1/2$ and got
$a(n+1) - (a(n) + a(n)^-1) = (-frac2sqrt8 + c) n^-3/2 + ldots$. So to get rid of the $n^-3/2$ term I want $c = frac2sqrt8$. Then look at the leading term for
$a(n) = sqrt2n + frac2sqrt8 ln(n) n^-1/2$ and continue in that vein...

From my answer here: Given $a_1=1, a_n+1=a_n+frac1a_n$, find $lim limits_ntoinftyfraca_nn$

Reposting here, as it is kind of lost in that thread and this thread is more suitable for it.

Note: I have no clue if a closed form exists, but here is an asymptotic estimate...

I think we can show that $$displaystyle a_n^2 sim 2n + dfraclog n2 - C$$ for some constant $displaystyle C gt 0$

By $displaystyle x_n sim y_n$ I mean $displaystyle lim_n to infty (x_n - y_n) = 0$

Consider $b_n = a_n^2 - 2n$

Then we have that $displaystyle b_n+1 = b_n + dfrac1b_n + 2n$

Notice that $b_0 gt 0$ and thus $displaystyle b_n gt 0$.

(Note that the other thread linked above starts with $a_1 = 1$ and not $a_0 = 1$.)

We can easily show that $b_n lt 2 + log n$, as

$b_n+1 - b_n = dfrac1b_n + 2n lt dfrac12n$

Adding up gives us the easy upper bound. Note, even though we can give tighter bounds, this is sufficient for our purposes.

Now we have that, for sufficiently large $displaystyle m,n$

$displaystyle b_m+1 - b_n = sum_k=n^m dfrac1b_k + 2k$

we have that

$displaystyle sum_k=n^m dfrac12k gt b_m+1 - b_n gt sum_k=n^m dfrac12k(1- dfracb_k2k)$

(Here we used $displaystyle dfrac11+x gt 1-x, 1 gt x gt 0$)

Now Since $b_k lt 2 + log k$, we have that

$displaystyle sum_k=n^m dfrac12k gt b_m+1 - b_n gt sum_k=n^m dfrac12k - sum_k=n^m dfrac2 + log k 4k^2$

Using the fact that $displaystyle H_m - H_n = log(dfracm+1n) + O(dfrac1n) + O(dfrac1n - dfrac1m)$, where $displaystyle H_n = sum_k=1^n dfrac1k$ is the $displaystyle n^th$ harmonic number.

We see that,

if $c_n = b_n - dfraclog n2$, then

$displaystyle O(dfrac1n -dfrac1m) + O(dfrac1n) gt c_m+1 - c_n gt O(dfrac1n -dfrac1m) + O(dfrac1n) -sum_k=n^m dfrac2 + log k 4k^2$

Now $displaystyle sum_k=1^infty dfrac2 + log kk^2$ is convergent and so by the Cauchy convergence criteria, we have that $displaystyle c_n$ is convergent.

Thus the sequence $displaystyle a_n^2 - 2n - dfraclog n2$ converges and hence, for some $displaystyle C$ we have that

$$displaystyle a_n^2 sim 2n + dfraclog n2 - C$$

or in other words

$$displaystyle a_n sim sqrt2n + dfraclog n2 - C$$

A quick (possibly incorrect) computer simulation seems to show a very slow convergence to $displaystyle C = 1.47812676429749dots$

Note: Didier suggested an alternate proof in the comments below, which might simpler.

Let us consider the functional formulation. Given $y(0)=1$, $y' = frac1y$ yields $ y(x) = sqrt2x+1$.

I am not saying $a(n)=y(n)$. Yet there is a link between the two approaches (finite difference).

Let $b_n=a_n^-1$, so that $b_n+1=(b_n+b_n^-1)^-1$ is the $n$th iterate of $$f(x)=(x+x^-1)^-1=x-x^3+x^5-x^7+cdots$$
The correct asymptotics have been given elsewhere, but this answer is to point out that there is a general method applicable to any function $f$ analytic at 0 with a series expansion $x+a_kx^k+a_k+1x^k+1+cdots$ with $k>1$ and $a_k<0$: I learned about it in Chapter 8 of de Bruijn's Asymptotic Methods in Analysis, where it is shown that if the starting value is small enough, the $n$th iterate of $f$ is asymptotically equivalent to $((1-k)a_kn)^-1/(k-1)$. In our case, $k=3$ and $a_k=-1$ so $b_n sim (2n)^-1/2$ and $a_n sim (2n)^1/2$.

A sketch of the proof in the case at hand is as follows. First, show by induction that if $u_n$ is any sequence tending to zero and such that
$$u_n+1=u_n-u_n^2+O(u_n^3)$$
then $u_n = n^-1 + O(n^-2log n)$. This is the tricky part of the proof.

Next make a substitution $z_n= 2b_n^2$. We have
$$ b_n+1=b_n(1-b_n^2+b_n^4-cdots)$$
so
$$z_n+1=z_nleft(1-frac12z_n + frac14z_n^2+ cdotsright)^2 = z_n-z_n^2 + frac34z_n^3 + cdots$$
from which we get $z_n = n^-1+O(n^-2log n)$ and the asymptotics for $b_n$ follow.

Incidentally, the sequence $b_n$ is what you get by using Newton's method to find a root of $xe^-x^-2/2$. But I couldn't find a way to make that shed any light on the asymptotics...