Explicit formula for the recurrence relation $x_n = a x_{n-1} + b (n-1) x_{n-2} +c$ Let $a,b,c$ be non-zero constants. How can we find an explicit formula for a sequence obeying a recurrence relation of the following type? $$ x_n = a x_{n-1} + b (n-1) x_{n-2} +c , \text{ for } n \geq 2$$ The initial conditions $x_0$ and $x_1$ are given, and can be assumed to be non-zero.
 A: An approach providing formulas, although not always very explicit ones, starts by converting the whole sequence $(x_n)$ into a generating function. The first idea that comes to mind to apply this technique might be to consider the simple generating function $$\sum_{n=0}^\infty x_nt^n$$ but, due to the $(n-1)x_{n-2}$ term in the recursion, one can predict that $x_n$ grows roughly as a factorial hence the radius of convergence of this series is $0$. 

The second idea that comes to mind might be to consider the exponential generating function $$X(t)=\sum_{n=0}^\infty x_n\frac{t^n}{n!}$$
Then the recursion yields $$X(t)=x_0+x_1t+\sum_{n=0}^\infty(a x_{n+1} + b (n+1) x_n +c)\frac{t^{n+2}}{(n+2)!}$$ that is, $$X(t)=x_0+x_1t+
a(X_1(t)-x_0t)+bX_2(t)+c\sum_{n=2}^\infty\frac{t^n}{n!}$$ with $$X_1(t)=\sum_{n=0}^\infty x_n\frac{t^{n+1}}{(n+1)!}\qquad
X_2(t)=\sum_{n=0}^\infty x_n(n+1)\frac{t^{n+2}}{(n+2)!}$$ One sees that $$ \sum_{n=2}^\infty\frac{t^n}{n!}=e^t-1-t$$
and that
$$
X_2(t)=\sum_{n=0}^\infty x_n(n+2-1)\frac{t^{n+2}}{(n+2)!}=tX_1(t)-Y(t)$$ with 

$$Y(t)=\sum_{n=0}^\infty x_n\frac{t^{n+2}}{(n+2)!}
$$

These definitions readily imply that
$$Y'(t)=X_1(t)\qquad X'_1(t)=X(t)$$ 
hence all these steps indicate that one should have considered from the start the exponential generating function $Y(t)$ rather than $X(t)$ since, translating everything in terms of $Y$, one gets 

$$Y''(t)-(a+bt)Y'(t)+bY(t)=Z(t)$$
  with $$Z(t)=x_0+x_1t-ax_0t+c(e^t-1-t)$$


The rest is standard. In the general case, the solutions of the homogenous equation $$U''(t)-(a+bt)U'(t)+bU(t)=0$$ span a two dimensional vector space. Pick a basis $\{U(t),V(t)\}$ of this vector space, then every solution $Y(t)$ of the complete differential equation is $$Y(t)=A(t)U(t)+B(t)V(t)$$ where $(A(t),B(t))$ solves the system $$A'(t)U(t)+B'(t)V(t)=0\qquad A'(t)U'(t)+B'(t)V'(t)=Z(t)$$ Thus, $$A'(t)=V(t)Z(t)W(t)\qquad B'(t)=-U(t)Z(t)W(t)$$ where $$W(t)=\frac1{U'(t)V(t)-V'(t)U(t)}$$ which yields $$Y(t)=U(t)\left(A_0+\int_0^tV(s)Z(s)W(s)ds\right)-V(t)\left(B_0+\int_0^tU(s)Z(s)W(s)ds\right)$$
where the constants $(A_0,B_0)$ are determined by the initial conditions $Y(0)=Y'(0)=0$, which read
$$U(0)A_0-V(0)B_0=0=U'(0)A_0-V'(0)B_0$$ hence, since $(U,V)$ is linearly independent, $A_0=B_0=0$, and, finally,

$$Y(t)=\int_0^t(U(t)V(s)-V(t)U(s))W(s)\,Z(s)\,ds$$

From there, depending on the exact form of the function $U$, $V$ and $W$, one can more or less easily deduce $Y(t)$, hence finally, every $x_n$ in terms of $(x_0,x_1)$.
