Attempting to solve this via generating functions. I decided to use exponential generating functions. Consequently, by shifting and index multiplying rules, we obtain
$$f(x)=xf(x)+ \int xf(x) dx$$
$$f'(x)=xf'(x)+f(x)+xf(x)$$
Solving this (initial condition $f(0)=1$) yields
$$f(x)= (e^{-x})/(1-x)^2 $$
So I used the following estimation to extract coefficients.
$$[z^n]\frac{f(z)}{(1-z)^\alpha}\sim f(1)\binom{n+\alpha-1}{n}\sim\frac{f(1)}{\Gamma(\alpha)}n^{\alpha-1}$$
That's where my issue is. If we use the first approximation, we get $(n+1)/e$ and if we use the second one then we get $n/e$.
To obtain the coefficients of the original sequence, I multiply by $n!$, since I used an EGF.
We obtain finally $(n+1)*n!/e$ or $n(n!)/e$.
From https://oeis.org/A000255, however, I was able to derive $(n+2)*n!/e$. This approximation agrees way more closely to the actual sequence than the above two (The error becomes less than $10^{-2}$ by the 8th term).
The error for the first two estimations though, only increases. Yes, the percent error itself decreases, because the difference between the approximations is multiplying by n or n+1 or n+2, which becomes less significant at greater values.
But I still feel that something's wrong. Are the first two estimations acceptable asymptotic versions produced by the formula, or did I mess up one of the many steps involved.
Cheers
The equation for the egf should be $$(x - 1) F''(x) + (x + 2) F'(x) + F(x) = 0,$$ with the general solution $$F(x) = \frac {C_1 e^{-x} + C_2 (x - 2) } {(x - 1)^2}.$$ Your solution corresponds to $f(0) = f(1) = 1$. Then the coefficient $[x^n] F(x)$ is asymptotically equivalent to $n/e$ (which is also a.e. to $(n + C)/e$ for any $C$).
$C = 2$ gives a better approximation because $2/e$ happens to be the next term in the asymptotic expansion of $[x^n] F(x)$. Since $$F(x) - \sum_{n \geq 0} \frac {n + 2} e x^n$$ extends to an entire function, we have $$[x^n] F(x) = \frac n e + \frac 2 e + O(\epsilon^n), \quad n \to \infty$$ for any $\epsilon > 0$.