How to Solve Advanced Recurrence Relations For some research I’m doing, I had to solve recurrences of the form: $a(n)=n \cdot a(n-1)-(n-3) \cdot a(n-2)$. 
WolframAlpha gave the solution:
$$ a(n) = 4c_2 \Gamma (n-1) ( \Gamma (n-1) F(1,n-2;n,n;1) -(c_1-4c_2(e-3))(n-2) )$$
Where F is the regularized generalized hypergeometric function.
I was taught a little bit about creating closed forms from recurrence relations, but have no clue to derive something such as this. It looks really similar to examples where you use exponential generating functions, but I can’t see exactly how it would work in this case. So my question is, how would a human get this solution, and where could I go to learn this and similar advanced techniques.
 A: 
Here are two references which might be useful when learning recurrence relations.
  
  
*
  
*Some interesting techniques are presented in Mathematics for the Analysis of Algorithms
  by D. E. Knuth and D. H. Greene. One of them is the so-called repertoire method. A real life application is presented in this post .
  
*Recurrences are also a prominent theme in Concrete Mathematics
  by R.L. Graham, D.E. Knuth and O. Patashnik.

Hint: Sometimes we are lucky and can find a solution also by elementary means. Assuming initial values $a_0, a_1$ are given, we derive from
\begin{align*}
a_n=na_{n-1}-(n-3)a_{n-2}\qquad n\geq 2
\end{align*}
successively
\begin{align*}
\color{blue}{1}a_2&=2a_1+a_0\\
a_3&=\color{blue}{3}a_2\\
a_4&=4a_3-a_2=\color{blue}{11}a_2\\
a_5&=5a_4-2a_3=\color{blue}{49}a_2\\
a_6&=6a_5-3a_4=\color{blue}{261}a_2\\
&\cdots\\
\end{align*}
Looking for the values $1,3,11,49,261$  in OEIS we find the sequence A001339 with representation
\begin{align*}
a_{n+2}=\sum_{k=0}^n\binom{n}{k}(k+1)!\qquad\qquad n\geq 0
\end{align*}
A: Hint.
Considering $S(x) = \sum_k a_k x^k$ we have
$$
S(x) -x\frac{d}{dx}\left(x S(x)\right)+x^4\frac{d}{dx}\left(\frac{S(x)}{x}\right) = 0
$$
and solving for $S(x)$ we obtain
$$
S(x) = C_0 (1-x)e^{-\frac 1x}
$$
A: Some general results that might be interesting.
As the OP said, generating functions do seem like a good method for this kind of recurrences. 
Let's consider a more general homogeneous second order linear recurrence with coefficients linear in $n$:
$$\left(A_2 n+B_2 \right) a_{n+2}+\left(A_1 n+B_1 \right) a_{n+1}+\left(A_0 n+B_0 \right) a_n =0$$
Now let's introduce some function defined by its Frobenius series:
$$f(x)=\sum_{n=0}^\infty b_n x^{n+s}$$
Which (we claim) obeys a second order ODE in the form:

$$(\alpha_2 x^2+\beta_2 x+\gamma_2) f''(x)+(\beta_1 x+\gamma_1) f'(x)+\gamma_0 f(x)=0$$

We have:
$$f'(x)=\sum_{n=-1}^\infty (n+1+s) b_{n+1} x^{n+s}$$
$$x f'(x)=\sum_{n=0}^\infty (n+s) b_n x^{n+s}$$
$$f''(x)=\sum_{n=-2}^\infty (n+1+s)(n+2+s) b_{n+2} x^{n+s}$$
$$x f''(x)=\sum_{n=-1}^\infty (n+s)(n+1+s) b_{n+1} x^{n+s}$$
$$x^2 f''(x)=\sum_{n=0}^\infty (n-1+s)(n+s) b_n x^{n+s}$$
Substituting the series, we obtain the following recurrence relation:
$$\alpha_2 (n-1+s)(n+s) b_n+\beta_2 (n+s)(n+1+s) b_{n+1} + $$ $$ + \gamma_2 (n+1+s)(n+2+s) b_{n+2} +\beta_1 (n+s) b_n+\gamma_1 (n+1+s) b_{n+1}+ $$ $$ +\gamma_0 b_n=0$$
Or:
$$\gamma_2 (n+s+1)(n+s+2) b_{n+2}+ $$ $$ +\left( \beta_2 n^2+(\beta_2 (2s+1)+\gamma_1)n+(\beta_2 s+\gamma_1) (s+1) \right) b_{n+1}+ $$ $$ +\left( \alpha_2 n^2+(\beta_1+(2s-1) \alpha_2)n+\alpha_2 s (s-1)+\beta_1 s+\gamma_0 \right) b_n=0 $$
To get rid of quadratic terms, we need to be able to factor out one of the brackets before the first term, let's pick $(n+s+1)$. Which means that the other two coefficients should have the form:
$$q(n+s+1)(n+p)=q(n^2+(p+s+1) n+p(s+1))$$
For the second term we have:
$$q=\beta_2 \\ p=s+\frac{\gamma_1}{\beta_2}$$
For the third term:
$$q=\alpha_2 \\ p=s+\frac{\beta_1}{\alpha_2}-2=\frac{\alpha_2 s (s-1)+\beta_1 s+\gamma_0}{\alpha_2 (s+1)}$$
We have obtained an equation which allows us to get rid of one parameter, for example, $\gamma_0$.
Now our recurrence relation is reduced to:

$$\gamma_2 (n+s+2) b_{n+2}+ \left(\beta_2n+\beta_2s+\gamma_1 \right) b_{n+1} +\left(\alpha_2 n+\alpha_2(s-2)+\beta_1 \right) b_n=0$$

We now have a linear second order recurrence with six parameters left, which (at least in principle) can be related to the original parameters $A_2,B_2,A_1,\dots$, and we can even set $a_n=b_n$.

In the example from Integral representation for the solution of a difference equation we have:
$$\gamma_2=1 \\ 
s=0 \\ 
\beta_2=\gamma_1=-2 \\
\alpha_2=1 \\
\beta_1=3 \\ \gamma_0=1$$
Which makes the generating function obey the following ODE:
$$(x-1)^2 f''(x)+(3 x-2) f'(x)+f(x)=0 \\ f(0)=1, \qquad f'(0)=0$$
Wolfram Alpha gives the exact solution:
$$f(x)=\frac{1}{1-x} \exp \left(- \frac{x}{1-x} \right)$$
From the explicit recurrence form given in the linked question, we also have:
$$f(x)= e \int_0^\infty e^{-(1-x) t} J_0 (2 \sqrt{t} ) dt$$
But this is correct and confirmed numerically.

So, the method allows us to obtain an ODE for a certain class of second order linear recurrence relations with variable coefficients linear in $n$.

Note: for the case in the OP the method doesn't seem to work as stated, some modifications need to be made.

Update: Another example.
Using hypergeometric equation, we can derive the following result. The sequence with explicit form:
$$a_n = \frac{(\alpha)_n}{(\beta)_n} {_2 F_1} (n+1, n+\alpha; n+\beta; p )$$
obeys the recurrence relation:
$$p(p-1) (n+2) a_{n+2}+ [(2p-1) n +(\alpha+2) p -\beta] a_{n+1} + (n+ \alpha) a_n=0 $$ $$ a_0 = {_2 F_1} (1, \alpha; \beta; p ) $$ $$ a_1 = \frac{\alpha}{ \beta}  {_2 F_1} (2, 1+\alpha; 1+\beta; p )$$
$p \neq 0, p \neq 1, \qquad \beta \notin \mathbb{Z}$.
For general initial conditions the solution also exists, but will be more complicated.
