Find an equivalent sequence as $n\to +\infty$ of $u_1>0, u_{n+1} = \frac{u_n}{n} + \frac{1}{n^2}$ Let $u_1>0$ be a real number. Let us consider $(u_n)_{n\geq 1}$ the sequence such as:
$$
\forall n \geq 1, u_{n+1} = \frac{u_n}{n} + \frac{1}{n^2}\quad (\star)
$$
Find an equivalent of $u_n$ as $n\to +\infty$.

So I found a way to show that $u_n \sim \frac{1}{n^2}$, but I'm quite unhappy with this method because I feel like I found it by chance without understanding anything (I did a lot of trials and found this)
My method: 
I showed by induction that $u_n \leq (u_1+1)$. Thus, $u_n\to 0$ considering $(\star)$.
Then, $nu_{n+1} = u_n + 1/n$. Thus (since $n+1 \sim n$), $nu_n \to 0$.
To end with, I have $n^2u_{n+1} = nu_n + 1$. Thus, $(n+1)^2 u_n \sim n^2u_n \to 1$.
How would you solve such a problem? Is there any more intuitive method that one may have done?
 A: 1. There is a heuristic argument which is useful for guessing the behavior of $u_n$: Rewrite the recurrence relation as
$$ u_{n+1} - u_n = \frac{u_n}{n} - u_n + \frac{1}{n^2}. $$
Its continuum analogue is the following differential equation:
$$ y' = \frac{y}{x} - y + \frac{1}{x^2}. $$
Using the standard method, this equation can be solved as:
$$ y(x) = x e^{-x} \int \frac{e^x}{x^3} \, \mathrm{d}x. $$
Then L'Hospital's Rule then tells that $y(x) \sim x^{-2}$ as $x \to \infty$. From this observation, we may as well expect that $u_n \sim n^{-2}$.
2. The above ansatz suggests that, in the recurrence relation for $u_{n+1}$, $\frac{1}{n^2}$ is the dominating term and $\frac{u_n}{n}$ is much smaller as $n\to\infty$. In particular, nesting this relation will produce an expansion with ever decreasing terms. This idea can be easily tested as follows:
Let $r_n = (n-1)u_n$. Then for $n \geq 2$,
$$ r_n = (n-1)u_{n} = \frac{1}{n-1} + \frac{r_{n-1}}{n-1}. $$
From this, we get
\begin{align*}
r_n
&= \frac{1}{n-1} + \frac{r_{n-1}}{n-1} \\
&= \frac{1}{n-1} + \frac{1}{(n-1)(n-2)} + \frac{r_{n-2}}{(n-1)(n-2)} \\
&= \frac{1}{n-1} + \frac{1}{(n-1)(n-2)} + \frac{1}{(n-1)(n-2)(n-3)} + \frac{r_{n-3}}{(n-1)(n-2)(n-3)} \\
&\qquad\vdots\\
&= \sum_{k=1}^{n-2} \frac{1}{(n-1)\cdots(n-k)} + \frac{r_2}{(n-1)!}
\end{align*}
Using this, it is not hard to conclude that $(n-1)^2 u_n \to 1$ as $n\to\infty$, and in fact, we can extract an asymptotic expansion of $u_n$ up to any prescribed order $\mathcal{O}(n^{-M})$.
A: This is similar to Sangchul Lee's approach. 
If we multiply $n!$ both sides of the recurrence, we obtain
$$
n!u_{n+1}=(n-1)!u_n+\frac{n!}{n^2}.
$$
Applying the above repeatedly, 
$$
u_{n+1}=\frac1{n!} \sum_{k=1}^n \frac{k!}{k^2} + \frac{u_1}{n!}.
$$
