Is there a general method to find the asymptotic order for this sequence? Given $$a_{n+1}=a_n+\frac{n}{a_1+\dots+a_n},\qquad a_1>0$$ The answer is $$\lim_{n\to\infty} a_n\sim\sqrt{3}\cdot\sqrt{n}-\frac{\sqrt{3}}{4}\cdot\frac{1}{\sqrt{n}}$$
It is easy to show this sequence is increasing, and is divergent, because if assume the opposite $A=\lim_{n\to\infty} a_n$, we will get $A=A+\frac{1}A$, which gives contradictions. We also have
$$a_{n+1}\ge a_n +\frac{1}{a_n}\ge2$$
From here I can show:
$$a_1+\frac{n}{a_n}\le a_{n+1}\le a_1+\frac{n}{a_1}$$
Update.(1)
This equation can be also written as:
$$\begin{align}
a_n&=\sum_{k=1}^n a_k-\sum_{k=1}^{n-1} a_k\\
\\
a_n&=\frac{n}{a_{n+1}-a_n}-\frac{n-1}{a_{n}-a_{n-1}}\\
\\
a_n(a_{n+1}-a_n)(a_{n}-a_{n-1})&=n(a_{n}-a_{n-1})-(n-1)(a_{n+1}-a_n)\tag{1}
\end{align}$$
If assume $a_n=c\cdot n^p$
$$a_{n+1}-a_n\sim cp\cdot n^{p-1},~~~a_{n}-a_{n-1}\sim cp\cdot n^{p-1}$$
Plug into $(1)$ and only keep the leading order:
$$c^2p\cdot n^{2p-1}=1$$
So we get $p=\frac{1}2$ and $c=\sqrt{2}$
Why the leading order coefficient is $\sqrt{2}$, not $\sqrt{3}$?
The bottom line is, if I take the answer as template, let $$a_n=c\sqrt{n}+t\frac{1}{\sqrt{n}}+o(\frac{1}{\sqrt{n}})$$
Pretend we don't know coefficients $c$ and $t$. Now we want to compute $c$ and $t$. To the leading order approximation, We have
$$a_{n+1}-a_n= c\cdot \frac{1}{2\sqrt{n}}+O(\frac{1}{n^{3/2}})$$
Next, plug into Eq.$(1)$ and we can solve for $c=\sqrt{2}$. Why does this give a contradiction?
Update.(2)
Thank you to @Sangchul Lee , @Somos and @Youem
I put the computation part in the answer box below, and it works for asymptotic approximation at arbitrary order.
 A: An alternative method depends on the sum sequence and a functional equation.
Define the sequence
$$ b_0 = 0,\;\; b_1 = a_1,\;\; b_2 = 2a_1 \!+\! 1/a_1,\; \text{ and }\\
b_n = 2 \!+\! (1 \!-\! b_{n-2}b_{n-3} \!-\! 2b_{n-2}^2)/b_{n-1}\; \text{ for }\; n>2. \tag{1}$$
This sequence satisfies the equation
$$ 0 = 1 + b_nb_{n+1} - 2b_{n+1}^2 + 2b_{n+2}^2 - b_{n+2}b_{n+3},
\quad \text{ for }\quad n\ge0. \tag{2}$$
Define the sequence $\,a_n := b_n-b_{n-1}.\,$ Then elementary algebra implies that
$$ a_{n+1} = a_n + \frac{n}{b_n} = a_n+\frac{n}{a_1+\dots+a_n},
\quad \text{ for }\quad n>0. \tag{3} $$
Assume the Ansatz
$$b_n=c\, n^p(1 + t_1x + t_2x^2 + t_3x^3 + \dots) \tag{4}$$ for some $\,p>0\,$ where $\,x := 1/n.\,$ Some algebra shows that
only $\,p=3/2\,$ works and moreover, for this value of $\,p,\,$
substituting the Ansatz into equation $(2)$ with $\,n+k\,$ replaced
with $\,x/(1+kx)\,$ gives
$\, 0 = 1-\frac34 c^2 + O(x).\,$ This implies
that $\,c=\frac2{\sqrt{3}}\,$ and more terms of the series expansion
for $\,b_n\,$ in terms of $\,x\,$ are found by solving for the
$\,t_k\,$ giving
$$ b_n = \frac2{\sqrt{3}}x^{-\frac32}\cdot\\
\left(1 \!+\! \frac{4c_0+3}4x \!+\!
\frac{8c_0^2 \!+\! 12c_0 \!+\! 3}{48}x^2 \!-\!
\frac{8c_0^3 \!+\! 18c_0^2 \!+\! 9c_0}{432}x^3 \!+\!
O(x^4)\right) \tag{5}$$
where $\,c_0\,$ is a constant depending on $\,a_1.\,$
The series expansion for $\,a_n\,$ is hence
$$ a_n = \sqrt{3}x^{-\frac12}\cdot\\
\left(1 + \frac{c_0}3 x - \frac{c_0^2}{18}x^2 + \frac{c_0^3}{54}x^3
- \frac{1520c_0^4 -81}{196992}x^4 +O(x^5)\right). \tag{6}$$

NOTE: The above analysis was purely formal. That is, I did not use
actual computation of $\,a_n\,$ for any particular values of $\,n\,$
and $\,a_1.\,$ However, I have now done the calculations of $\,a_n\,$
where $\,n=1\,$ up to $\,n=2^{15}\,$ and for $\,a_1=1\,$ up to $\,a_1=10.\,$
Assuming $\,a_n = \sqrt{3n}(1+c_n/n)\,$ where $\,c_n\,$ depends on
$\,a_1,\,$ then plotting $\,c_n\,$ with $\,y=c_n\,$ versus $\,x=\log(n)\,$
seems to show that $\,y\,$ oscillates very roughly between
$\,-.14\,a_1^2\,$ and $\,.11\,a_1^2\,$ with a period of roughly $9.$
Thus, $\,c_n\,$ does not have a limiting value and it
does depend on $\,a_1.\,$ This is all based on limited
numerical data and is not even close to a proof, still it
is suggestive of the true situation.
A: 
$$a_{n+1}=a_n+\frac{n}{\sum_{k=1}^n a_k}$$

This equation can be written as:
$$\sum_{k=1}^n a_k=\frac{n}{a_{n+1}-a_n}$$
Further, we can find $a_n$ by taking subtraction:
$$\begin{align}
a_n&=\sum_{k=1}^n a_k-\sum_{k=1}^{n-1} a_k\\
\\
a_n&=\frac{n}{a_{n+1}-a_n}-\frac{n-1}{a_{n}-a_{n-1}}\\
\\
a_n(a_{n+1}-a_n)(a_{n}-a_{n-1})&=n(a_{n}-a_{n-1})-(n-1)(a_{n+1}-a_n)\tag{1}\\
\\
&=(a_{n+1}-a_n)-n(a_{n+1}-2a_n+a_{n-1})
\end{align}$$
First determine the leading order. Let $a_n=c_0\cdot n^p$
$$\begin{align}
a_{n+1}&=c_0(n+1)^p=c_0 \left(n^p+pn^{p-1}+\frac{1}{2}p(p-1)n^{p-2}+...\right)\tag{2}\\
\\
a_{n-1}&=c_0(n-1)^p=c_0 \left(n^p-pn^{p-1}+\frac{1}{2}p(p-1)n^{p-2}+...\right)\tag{3}
\end{align}$$
Plug into $(1)$ and only keep the leading order on both sides:
$$c_0n^p\cdot c_0pn^{p-1}\cdot c_0pn^{p-1}=(2 c_0 p - c_0 p^2)n^{p-1}$$
So we get:
$$3p-2=p-1,~~~~~~c_0^3 p^2=(2 c_0 p - c_0 p^2)$$
$$\Rightarrow p=\frac{1}2,~~~c_0=\sqrt{3}$$
Next, find the next order, so we set the solution Ansatz as:
$$a_n=\sqrt{n}\left(c_0+c_1\cdot\frac{1}n \right),~~~~~c_0=\sqrt{3}$$
Plug into $(1)$ and re-do the series expansion to update $(2)$ and $(3)$. Now, comparing the coefficients for $\frac{1}{n^{3/2}}$
$$-\frac{3}4c_1\frac{1}{n^{3/2}}=\left(-\frac{\sqrt{3}}8 - \frac{5}4c_1\right)\frac{1}{n^{3/2}}$$
$$\Rightarrow c_1=-\frac{\sqrt{3}}{4}$$
Next, set the solution Ansatz as:
$$a_n=\sqrt{n}\left(c_0+c_1\cdot\frac{1}n+c_2\cdot\frac{1}{n^2} \right),~~~~~c_0=\sqrt{3},~~c_1=-\frac{\sqrt{3}}{4}$$
Now, comparing the coefficients for $\frac{1}{n^{5/2}}$
$$\frac{3}{32}(\sqrt{3}-40c_2)\frac{1}{n^{5/2}}=\frac{3}{64}\left(\sqrt{3} - 112c_2\right)\frac{1}{n^{5/2}}$$
$$\Rightarrow c_2=-\frac{\sqrt{3}}{32}$$
Keep going on this method, we get:
$$a_n=\sqrt{n}\left(\sqrt{3}-\frac{\sqrt{3}}{4}\cdot\frac{1}n-\frac{\sqrt{3}}{32}\cdot\frac{1}{n^2}-\frac{\sqrt{3}}{128}\cdot\frac{1}{n^3}-\frac{79\sqrt{3}}{38912}\cdot\frac{1}{n^4}+\cdots \right)$$
