<s>[INST] Recurrence Relation How do I solve:
$k(k+1)a_{k}=2(\lambda k-1)a_{k-1}+(a-\lambda^2)a_{k-2}$
where $\lambda$ and $a$ are constants, and similar other recurrence relations?
 [/INST] Alright, let us substitute for $k-2$ . then the recurrence relation reads:
\begin{equation}
(k+2)(k+3) a_{k+2} - 2 \left( \lambda (k+2)-1 \right) a_{k+1} - \left( a+\lambda^2 \right) a_k = 0
\end{equation}
for $k=0,1,2\cdots$. Here for the sake of simplicity we will make some additional assumption that $a_1 = (\lambda-1) a_0$. We will present the solution in this case and then we will show how to wave this assumption.  Let us get it over and done with now.
We define the Z-transform as $A[z]:= \sum\limits_{n=0}^\infty a_n z^{-n}$. Then simple algebra shows that the quantity in question satisfies the following ODE:
\begin{equation}
\frac{d^2 A[z]}{d z^2} + \frac{2\lambda}{z^2} \frac{d A[z]}{d z} - \frac{a+\lambda^2-2 z+2 \lambda z}{z^4} \cdot A[z] = 0
\end{equation}
Now we use standard techniques to eliminate the coefficient at the first derivative. We have $A[z] = \exp(-1/2 \int (2 \lambda/z^2) \quad dz) \cdot v(z) = \exp(\lambda/z) \cdot v(z)$ and the function $v(z)$ satisfies the following ODE:
\begin{equation}
\frac{d^2 v(z)}{d z^2} - \frac{a+\lambda^2-2 z}{z^4} v(z) = 0
\end{equation}
It is not hard to see that the ODE above can be obtained from the confluent hypergeometric ODE by changing both the abscissa and the ordinate accordingly. The final result is given below:
\begin{eqnarray}
a_n &:=& \left. \frac{1}{n!} \frac{d^n }{d x^n} 
F_{1,1} \left[1- \frac{1}{\sqrt{a+2 \lambda^2}}, 2, 2 \sqrt{a+2 \lambda^2} x \right] \cdot \exp\left( (\lambda-\sqrt{a+2 \lambda^2}) x\right)
\right|_{x=0}\\
&=&
\frac{2^n \left(a+2\lambda^2\right)^{n/2}  \Gamma \left(n-\frac{1}{\sqrt{2\lambda^2+a}}+1\right)}{\Gamma (n+1) \Gamma (n+2) \Gamma
   \left(1-\frac{1}{\sqrt{2\lambda^2+a}}\right)}
\, _2F_1\left(-n-1,-n;\frac{1}{\sqrt{2\lambda^2+a}}-n;\frac{1}{2}-\frac{\lambda}{2 \sqrt{2\lambda^2+a}}\right)
\end{eqnarray}
for $n=0,1,2,\cdots$.
The Mathematica code below demonstrates that the results are correct:
In[862]:= {a, l, z} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50];
aa = Table[((2^n) ((a + 2 l^2)^(n/2))
       Gamma[1 - 1/Sqrt[a + 2 l^2] + n] )/(
     Gamma[1 - 1/Sqrt[a + 2 l^2]] Gamma[1 + n] Gamma[2 + n])
      Hypergeometric2F1[-1 - n, -n, 1/Sqrt[a + 2 l^2] - n, 
      1/2 - l/(2 Sqrt[a + 2 l^2])], {n, 0, 10}] // Simplify;
(aa[[2]] - (l - 1) aa[[1]]) // Simplify
Table[(k + 2) (k + 3) aa[[k + 3]] - 
   2 (l (k + 2) - 1) aa[[k + 2]] - (a + l^2) aa[[k + 1]], {k, 0, 
   Length[aa] - 3}] // Simplify

Out[864]= 0.*10^-50

Out[865]= {0.*10^-49, 0.*10^-50, 0.*10^-50, 0.*10^-50, 0.*10^-50, 
 0.*10^-51, 0.*10^-51, 0.*10^-52, 0.*10^-52}

</s>