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 good first approach for many such problems is to try ansätze similar to what would solve a corresponding differential equation. $y' = (ax+b)y$ has a solution of the form $y = c_1^{(x^2)}c_2^x \equiv e^{(c_1x^2+c_2x)}$. Adding the constant $+c$ gives a complicated diffeq with erfc in there. Since it's a second-order recurrence it's probably also good to try a sum of two such terms. – Alex Meiburg Jul 28 '17 at 7:43
• I found that oeis.org/A000932 and oeis.org/A059480 have the explicit solutions (for initial conditions $(0,1)$ and $(1,0)$ which is sufficient to construct any) for $a=b=1$ in terms of Hypergeometric functions. The fact that the characteristic polynomial there "isn't anything in particular" leads to me to believe that it could be generalized to arbitrary $a$ and $b$. – Alex Meiburg Jul 28 '17 at 17:58

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)$.