Consider a following recursion relation of degree two: \begin{equation} y_{n+1} = f_n \cdot y_n + g_n y_{n-1} \end{equation} for $n\ge 1$ subject to $y_1={\mathcal F}_1$ and $y_0={\mathcal F}_0$. By backwards iterating this recurrence we have shown that the solution reads: \begin{equation} y_{n+1} = (\prod\limits_{q=1}^n f_q) \sum\limits_{l=0}^{\lfloor \frac{n+1}{2} \rfloor} \sum\limits_{l+2 \le p_l < p_{l-1} < \cdots < p_1 \le n+2} \left(\prod\limits_{\eta=1}^l \frac{g_{p_\eta-\eta-1}}{f_{p_\eta-\eta-1}f_{p_\eta-\eta-2}}\right) \cdot (1_{l+2 < p_l} {\mathcal F}_1 + 1_{l+2 = p_l} {\mathcal F}_0) \end{equation} The formula above can be brought into closed form if we assume that: \begin{equation} \frac{g_n}{f_n (1_{n > 1} f_{n-1} + 1 \delta_{n,1})} = const(n)=\theta \end{equation} In that case we get: \begin{eqnarray} &&\frac{y_{n+1}}{(\prod\limits_{q=1}^n f_q)}= \\&& \left( ({\mathcal F}_0 -{\mathcal F}_1)\cdot \theta \, _2F_1\left(\frac{1}{2}-\frac{n}{2},1-\frac{n}{2};1-n;-4 \theta\right)+ {\mathcal F}_1 \, \cdot _2F_1\left(-\frac{n}{2}-\frac{1}{2},-\frac{n}{2};-n-1;-4 \theta\right) \right) \end{eqnarray} If we now use the algorithm Hyper from the $A=B$ book and then we get: \begin{equation} \frac{y_{n+1}}{(\prod\limits_{q=1}^n f_q)}= (-1)^{-n} \frac{ \left(({\mathcal F}_0 z_{+}+{\mathcal F}_1) \left(z_{-}^{n+1}-z_{+}^{n+1}\right)-{\mathcal F}_0 z_{+}^{n+1} (z_{-}-z_{+})\right)}{z_{-}-z_{+}} \end{equation} where $z_{\pm} := \frac{1}{2} \left(-1\pm \sqrt{1+4 \theta}\right)$.
All what was above did not use any assumptions on the sequence $f_n$. Now let us assume that $f_n:=(n+a)/(n+b)$ for some $a,b \in {\mathbb Q}$. Then we get: \begin{eqnarray} &&{\mathcal F}_1 + \sum\limits_{n=1}^\infty y_{n+1} z^n = \\ && \frac{({\mathcal F}_1 + {\mathcal F}_0 z_+) z_- F_{2,1}\left[\begin{array}{rr} 1 & 1+a \\ 1+b \end{array};-z_- z\right] - ({\mathcal F}_1+{\mathcal F}_0 z_-) z_+ F_{2,1} \left[ \begin{array}{rr} 1 & 1+a \\ 1+b \end{array};-z_+ z\right]}{z_--z_+} \end{eqnarray} whenever the series on the left hand side are convergent.
In[3821]:= th = RandomInteger[{1, 50}]; th = -th/200;
{a, b} = RandomInteger[{1, 20}, 2];
{F0, F1} = RandomInteger[{1, 20}, 2];
f[n_] := (n + a)/(n + b);
g[n_] := th f[n] If[n == 1, 1, f[n - 1]];
{zp, zm} = {1/2 (-1 + Sqrt[1 + 4 th]), 1/2 (-1 - Sqrt[1 + 4 th])};
M = 200;
y = Table[0, {M}]; y[[{1, 2}]] = {F0, F1};
For[n = 1, n <= M - 2, n++,
y[[n + 2]] = f[n] y[[n + 1]] + g[n] y[[n]];
];
z = RandomReal[{0, 1}, WorkingPrecision -> 50];
Total[Table[y[[n + 2]] z^n, {n, 1, M - 2}]]
1/(zm - zp) ((F1 + F0 zp) zm Hypergeometric2F1[1, 1 + a,
1 + b, -zm z] - (F1 + F0 zm) zp Hypergeometric2F1[1, 1 + a,
1 + b, -zp z]) - F1
Out[3831]= 788.290548731724000915746764646576579876823076494
Out[3832]= 788.290548731724000915808972672164379808631261417
Having said all that above my question would be in what other cases can we do the multiple sum in the second equation from the top and therefore obtain closed form expressions for our recursion ?
Let us now assume that: \begin{equation} \frac{g_n}{f_n f_{n-1}} = \Theta \cdot \binom{n}{\theta} \end{equation} where $\theta$ is a fixed positive integer.Then using the second equation from the top we get: \begin{eqnarray} (I) \quad\quad \frac{y_n}{(\prod\limits_{q=1}^n f_q)} = \left( ({\mathcal F}_0-{\mathcal F}_1) \cdot (\sum\limits_{l=1}^{\lfloor \frac{n+1}{2} \rfloor} \Theta^l {\mathcal S}^{(l+1,n)}_{l-1}(\theta)) + {\mathcal F}_1 \cdot (\sum\limits_{l=0}^{\lfloor \frac{n+1}{2} \rfloor} \Theta^l {\mathcal S}^{(l,n)}_{l}(\theta)) \right) \end{eqnarray} where \begin{eqnarray} {\mathcal S}^{(m,n)}_l(\theta) &:=& \sum\limits_{m+2 \le p_l < p_{l-1} < \cdots < p_1 \le n+2} \prod\limits_{\eta=1}^l \left(p_\eta-\eta-1\right)_{(\theta)} \\ &=& \sum\limits_{\lambda=0}^{\theta \cdot (l-1)} \binom{n+2-l}{l+\theta+\lambda} \binom{(l-1) \theta}{\lambda} \cdot \sum\limits_{j=0}^{(\theta-1)\cdot (l-1)} {\mathcal B}^{(l-1)}_j (\theta) \cdot \lambda_{(j)} \end{eqnarray} for $\theta > m-l$. Here the coefficients on the right hand side satisfy following recursion relations: \begin{eqnarray} {\mathcal B}_J^{(d-1)}(\theta) = \sum\limits_{j=0\vee (J-\theta)}^{J \wedge [(d-2)\cdot (\theta-1)]} {\mathcal B}_j^{(d-2)}(\theta) \cdot (-1)^{J-j} \binom{2 d-2+J+\theta}{\theta} \binom{\theta}{J-j} \cdot \frac{\left((d-2) \theta-j+1\right)^{(j)}}{\left((d-1) \theta-J+1\right)^{(J)}} \cdot \frac{\left(j-(d-1)(\theta-1)\right)^{(J-j)}}{\left(j+2 d-1 +\theta\right)^{(J-j)}} \end{eqnarray} subject to ${\mathcal B}^{(0)}_0=1$.
Now, three notes are in place.
Firstly the solution (I) to the original recursion relation consists of two modes the first one being proportional to the "initial velocity" (${\mathcal F}_0 - {\mathcal F}_1$) and the second one to the initial position (${\mathcal F}_1$) all this just like in the case of a order-two recursion relation with constant coefficients.
Secondly, this example is not solvable in the sense of the algorithm Hyper (maybe except for the case $\theta=1$--please check it) however our solutions are also "closed form" in the sense that the computational expense is proportional to $n^2$.
Thirdly, let us try to solve the last recursion relation for the coefficients ${\mathcal B}^{(d-1)}_J$. Let us fix $J$. Then the sum on the rhs runs over $j=0,\cdots, J$. Assume by induction that we know all coefficients ${\mathcal B}^{(d-1)}_j$ for $j=0,\cdots,J-1$. Then clearly the recursion is just a first order recursion with a source term on the rhs. As such that recursion can be always "solved" as a convolution of the relevant coefficient with the rhs. It is just not clear whether those convolutions can be expressed in closed form. We carried out this procedure for $J=0,1$ and we got: \begin{eqnarray} {\mathcal B}_0^{(d)} &=& \prod\limits_{\eta=0}^d \binom{2 \eta+\theta}{\theta} \\ {\mathcal B}_1^{(d)} &=& \left(\frac{\binom{d+\frac{\theta }{2}+\frac{1}{2}}{d+1}}{\binom{d+\frac{1}{2}}{-\frac{1}{2}}}-(d+1) \theta -1\right) \frac{(\theta -1) }{d (\theta -2) \theta} \cdot \prod\limits_{\eta=0}^d \binom{2 \eta+\theta}{\theta} \end{eqnarray} Now, when we go to $J=2$ the situation gets messy. The formula has the same structure as in the two cases listed above but the pre-factor now depends on harmonic numbers and as such is not a hyper-geometric function . Thus we conclude that is only for $\theta=1$ that the solution has a closed form. In all other cases we need to pre-calculate the ${\mathcal B}^{(d)}_J$ coefficients and use them in calculations.