How to solve this recurrence Relation - Varying Coefficient Sir,I have two questions related to this recurrence relation. It has been messing with me for long. Because of this I couldn't proceed my work for some time .This contains a polynomial term n+2 in sequence. I tried almost all known methods to me. I tried generating functions,general methods all etc . So it makes difficult for solving it. Would you provide the solution for it
$ R_{n}=\frac{1}{n} \{C_1 R_{n-1} +C_2 R_{n-2}\}, R_0 = A_0 ,R_1 = A_1  $
Here $C_1$  and $C_2 $ are constants 
Question 1 ::a) What is the solution of this in terms of n for $R_n $ b) What is the summation($ \sum_{n=0}^{n= \infty}  R_n$) of it, in-terms of n if possible provided $R_0,R_1,C_1,C_2, R_n  $ are numerical constants of singular matrices. Means normal case,numerical constants
Question 2 ::a) What is the solution of this in terms of n for $R_n$ and What is the summation($ \sum_{n=0}^{n= \infty}  R_n$) of it, in-terms of n if possible provided $ R_0,R_1,C_1,C_2 $ are   matrices of order $3\times 3 $?
 A: For now I'll treat only the case in which the terms of the recurrence are scalars. My strategy will be to introduce the generating function $R(x) = \sum_{n=0}^\infty$, simply it as best as possible, and then examine $R(1)=\sum_{n=0}^\infty$. Rewriting your recurrence relation slightly, we want to solve $$(n+2) R_{n+2} = C_1 R_{n+1}+C_2 R_{n}.$$ If we sum over $x^n$ on both sides we obtain the sums
\begin{align}
\sum_{n=0}\,(n+2) R_{n+2} x^n
&= \sum_{n=2}^\infty\, n R_n x^{n-2} \\
&= \frac{1}{x^2}\left[\left(\sum_{n=0}^\infty\,n R_n x^n\right)-R_1 x\right]\\
&=\frac{1}{x^2}\left[x \frac{d}{dx}\left(\sum_{n=0}^\infty R_n x^n\right)-R_1 x\right] \\
&= \frac{1}{x}\left[R'(x) - R_1 \right],
\end{align}
$$\sum_{n=0}^\infty R_{n+1} x^n =\frac{1}{x} \sum_{n=0}^\infty  R_{n+1}x^{n+1} =\frac{1}{x}\left[R(x)-R_0\right].$$ Therefore $$\frac{1}{x}\left(R'(x) - R_1 \right) = \frac{C_1}{x}\left[R(x)-R_0\right]+C_2 R(x),$$ which gives the inhomogeneous ODE $$R'(x)-(C_1 +C_2 x) R(x) = R_1-C_1 R_0.$$ By introducing an appropriate integrating factor we can rewrite this as $$\frac{d}{dx}\left(e^{-C_1 x-C_2 x^2/2}\right)=e^{-C_1 x- C_2 x^2/2}\left(R_1-C_1 R_0 \right)$$  and then integrate both sides to obtain the formal solution
\begin{align}
R(x)
&=e^{C_1 x+C_2 x^2/2}\left[R_0+\left(R_1-C_1 R_0\right)\int_0^x dy\, e^{-C_1 y- C_2 y^2/2} \right]\\
&=e^{C_1 x+C_2 x^2/2}\left[R_0+\left(R_1-C_1 R_0\right) e^{C_1^2/2 C_2} \sqrt{\frac{\pi }{2 C_2}} \left(\text{erf}\left(\frac{C_1+C_2 x}{\sqrt{2 C_2}}\right)-\text{erf}\left(\frac{C_1}{\sqrt{2 C_2}}\right)\right)\right] \end{align}
The summation is now easy to extract as
\begin{align}
R(1)
&=e^{C_1+C_2/2}\left[R_0+\left(R_1-C_1 R_0\right) e^{C_1^2/2 C_2} \sqrt{\frac{\pi }{2 C_2}} \left(\text{erf}\left(\frac{C_1+C_2}{\sqrt{2 C_2}}\right)-\text{erf}\left(\frac{C_1}{\sqrt{2 C_2}}\right)\right)\right] \end{align}
This represents the final answer in the scalar case (at least discounting possible algebra errors...)
What about generalizing to $R$ and $C$ being $N\times N$ matrices? I haven't as much expertise in this realm, but I think that everything remains the same up to the ODE (I see no matrix multiplication going on). The formal solution should be the same, except that now the integrating factors are matrix exponentials in $(C_1,C_2)$. So in principle 'all' one needs to do is compute these exponentials... (See Reyal's newer question for discussion of such.)
