Solution of second order difference equation with non constant coefficient What is the general approach to solve a second order difference equation with non constant coefficients? 
e.g $$X_{n+1} = 5nX_n + X_{n-1}$$
 A: Use generating functions. Define $g(z) = \sum_{n \ge 0} X_n z^n$, shift your recurrence, multiply by $z^n$ and add over $n \ge 0$, recognize resulting sums:
$\begin{align*}
  X_{n + 2}
    &= 5 (n + 1) X_{n + 1} + X_n \\
  \sum_{n \ge 0} X_{n + 2} z^n
    &= 5 \sum_{n \ge 0} (n + 1) X_{n + 1} z^n + \sum_{n \ge 0} X_n z^n \\
  \frac{g(z) - X_0 - X_1 z}{z^2}
    &= 5 g'(z) + g(z)
\end{align*}$
This is a linear differential equation of the first order for $g$, extract coefficients from its solution to get $X_n$
A: The following is motivated by S. Janson's paper "A divergent generating function that can be summed and analysed analytically" cited in "Formulas" for OEIS sequence A001040. The link given at the OEIS page does not seem to work, but the paper is available at the author's website here.
The modified Bessel function of the first kind, denoted $I_n(z)$, is known to satisfy the recurrence relations
$$I_{n-1}(z)-I_{n+1}(z)=\frac{2n}{z} I_{n}(z)$$
(see for instance DLMF 10.29.1, with notation defined in DLMF 10.25). If we choose $z=-2/5$ and rearrange slightly, we obtain
$$I_{n+1}(-2/5)=5n I_{n}(-2/5)+I_{n-1}(2/5).$$
Hence the recurrence relation for $I_n(-2/5)$ coincides with that of $X_n$, so $X_n=I_n(-2/5)$ is one solution to the recurrence relation. Similarly the modified Bessel function of the second kind, denoted $K_n(z)$, satisfies the recurrence relation
$$K_{n-1}(z)-K_{n+1}(z)=-\frac{2n}{z} K_{n}(z).$$
In the case of $z=2/5$ this rearranges to
$$K_{n+1}(2/5)=5n K_{n}(2/5)+K_{n-1}(2/5)$$ and we conclude that $X_n=K_{n+1}(2/5)$ is also a solution to the recurrence relation. Since these two solutions are linearly independent, the general solution to $X_{n+1}=5n X_n+X_{n-1}$ is of the form $$X_n=a I_n(-2/5)+b K_n(2/5).$$We can in fact replace $5n$ with any linear function $p+qn$ by shifting $n\to n+p/q$ in the Bessel recurrence relations and choosing $z=-2/$, yielding $$X_n = a I_{n+p/q}(-2/q)+b K_{n+p/q}(-2/q)$$ as the general solution of $X_{n+1}=(p+qn)X_{n}+X_{n-1}$.
It is thus possible to provide solutions to non-constant recurrence relations of this kind. But the reader may rightly complain about the artificiality of this solution: How would we have known to consider modified Bessel functions in the first place? The details are substantial and so I'll largely leave them to the paper of S. Janson. The basic idea is to consider the (formal) generating function $x(t)=\sum_{n=0}^\infty X_n t^n$. By multiplying both sides of the recurrence relation by $t^{n-1}$ and summing from $n=1$ to $\infty$, we get an ODE in $x(t)$. In principle, then, the general form of $X_n$ corresponds to coefficients in the general solution to this ODE.
The tricky point here is that, for almost initial conditions, the sequence $\{X_n\}$ will diverge rapidly---so rapidly, in fact, that the function $x(t)$ will have zero radius of convergence. (More precisely, the large-order asymptotics of $K_\nu(z)$ yields $X_n\sim b K_{n+p/q}(2/a)\sim \frac12 b(n-1)!a^{n+p/q}n^{p/q}$.) To handle $x(t)$ sensibly therefore more subtle tools like Borel summation, but that's about where my patience runs out and I'll pursue it no further.
