The asymptotics of a recursion I am trying to assess the aymptotic behavior for large $n$ of the sequence $(a_n)_{n=1}^\infty$ defined by the recursion
\begin{align}
a_{n} &= (n+1) a_{n-1} - n a_{n-2}-1, \ \forall n\ge3 \\ 
a_2 &= 2a_1-3.
\end{align}
This recursion originates from this question on stats.stackexchange.com.
Define the generating function $f(x):=\sum_{n=0}^\infty a_nx^n,$. $g(x):=x^2f(x)$ satisfies the ODE
$$g'-\Big(\frac1{x^2}+\frac1x+\frac1{1-x}\Big)g=\Big(\frac x{1-x}\Big)^2-\frac{a_0}{1-x}+(2a_0-a_1)\frac x{1-x}$$
Multiply both sides by $e^h=\frac{1-x}xe^{\frac1x}$ where $h(x):=\frac1x+\ln\frac{1-x}x$ and get
$$\frac{d}{dx}\Big(\frac{1-x}xe^{\frac1x}g\Big)=e^{\frac1x}\Big(\frac x{1-x}-\frac{a_0}x+(2a_0-a_1)\Big).$$
Now $e^{\frac1x}$ seems to prevent me from solving the ODE in a "closed" form. 
I read somewhere the complex analysis may help to derive the large $n$ asymptotics of $a_n$. How does one procede?
 A: Partial answer to an interesting problem.  However, without $a_0$ and $a_1$ specified, one can't even determine the sign of the asymptotics.  By dropping the -1 on the right-hand side of the equation, the recursion can be solved (in Mathematica) in terms of 'subfactorials.'  That is
$$a_n = \kappa_1+ \kappa_2 \sum_{k=0}^n k! $$ 
where the $\kappa$ are related somehow to $a_0$ and $a_1.$  Naturally the second term will dominate.  Choosing $a_0=0$ and using numerical and symbolic computations I was able to determine
$$ a_n \sim (a_1+2-e)\sum_{k=0}^n k! \quad ,\quad (a_0 = 0)$$
To make it absolutely clear, this step lacked rigor. The sum may be approximated by
$ \sum_{k=0}^n k! \sim n! (n-1)/n .$  As a check, for $a_1=2,$ by the time $n=10,$ the relative error is about 1 part in a million.  What would be an interesting study is what happens to the asymptotics when $a_1 = e-2;$ it appears to be slowly growing, like $\cal{o}(n^{1/4}).$  We therefore have a region where, even with one constant simply specified, a non-uniformity appears.  With two free constants, this problem might be a real bear to get a handle on.
Update: The second formula is a proposed asymptotic solution for
$a_n = (n+1)a_{n-1}-n \ a_n -1$ for $a_0=0$ and away from the region $a_1=e-2.$  The first asymptotic form was deduced by using Mathematica's RSolve command.  To check, and to generate numerical data for my guess, I did the following:
$$\text{ Clear[a]; a[0]=0; a[1]=1; a[n_]:=a[n]=(n+1)a[n-1]-n*a[n-2]-1;}$$ 
$$\text{tb=Table[a[n],{n,0,30}]/Table[ Sum[k!,{k,0,n}],{n,0,30}]}$$
The output is 
$$ \text{N[tb]={0, .5, .5, .4, .323529, ... 0.281718}} $$
One of the last numbers looked a lot like $e$ minus a rational, so I used 
$$ \text{FindIntegerNullVector[{1,E,0.281718...}]}$$
to confirm my suspicion.  Doing this for many $a_1$ led to the conjecture; it is not a proof I don't have enough motivation to go further.
A: Define $b_n:=a_n-a_{n-1}$. The original recursion turns into
$$b_n=nb_{n-1}-1$$
which is solvable with a generating function. 
