Finding the closed form of two intertwined recurrence relations I would like to find the closed form (if there exists one) of the two recurrence relations
\begin{align*}
a_n&=(\gamma-n+1)a_{n-1}+\omega b_{n-1}\\
b_n&=(\gamma-n+1)b_{n-1}-\omega a_{n-1},
\end{align*}
with the initial conditions $a_0=1$ and $b_0=0$, where $\gamma\in(0,1)$ and $\omega>0$. The first few terms of the sequences are
\begin{align*}
a_0&=1 & b_0&=0\\
a_1&=\gamma & b_1&=-\omega\\
a_2&=\gamma(\gamma-1)-\omega^2 & b_2&=-\omega(\gamma-1)-\omega\gamma\\
\end{align*}
I've tried using summations to find the solution. Let $A(x)=\sum_{n=0}^\infty a_nx^n$ and $B(x)=\sum_{n=0}^\infty b_nx^n$. Then
\begin{align*}
\sum_{n=1}^\infty a_nx^n&=\sum_{n=1}^\infty(\gamma-n+1)a_{n-1}x^n+\sum_{n=1}^\infty\omega b_{n-1}x^n\\
&=\gamma\sum_{n=1}^\infty a_{n-1}x^n-\sum_{n=1}^\infty(n-1)a_{n-1}x^n+\omega\sum_{n=1}^\infty b_{n-1}x^n\\
&=\gamma x\sum_{n=0}^\infty a_nx^n-x^2\sum_{n=0}^\infty na_nx^{n-1}+\omega\sum_{n=0}^\infty b_nx^n\\
&=\gamma xA(x)-x^2A'(x)+\omega B(x).
\end{align*}
which implies
$$
A(x)-1=x\gamma A(x)-x^2A'(x)+\omega xB(x).
$$
Similarly, we have
$$B(x)=x\gamma B(x)-x^2B'(x)-\omega xA(x).$$
We can then use substitution to get an extremely nasty second-order differential equation. Is there any way to go from here, or is there another method that is known on how to solve recurrence relations like these?
 A: Notice that 
$$\omega b_{n-1}=a_n-(\gamma-n+1)a_{n-1}.$$
Substituting into the second equation yields the recurrence relation
\begin{align}a_n=(2\gamma-2n+3)a_{n-1}-((\gamma-n+2)^2+\omega^2)a_{n-2}, \tag{$\ast$}\end{align}
with $a_0=1$ and $a_1=\gamma$.
Using generating function $A(x)=\sum_{n=0}^\infty a_nx^n$, we get the differential equation
$$\small{ x^4A''(x)+(x^3-2\gamma x^2+2x)A'(x)+((\gamma^2+\omega^2)x^2-(2\gamma+1)x+1)A(x)+2(\gamma+1)x-(1+\gamma)=0}.$$
This is beyond my abilities to solve. However, Wolfram Alpha solves ($\ast$), giving
$$a_n=\frac12(-1)^n\left(\frac{\Gamma(n-(\gamma+i\omega))}{\Gamma(-(\gamma+i\omega))}+\frac{\Gamma(n-(\gamma-i\omega))}{\Gamma(-(\gamma-i\omega))}\right).$$
Since $n\in\mathbb{N}$, we get
\begin{align*}
a_n&=\frac12((\gamma+i\omega)_n+(\gamma-i\omega)_n)\\
&=\frac12\sum_{k=0}^n\left[s(n,k)(\gamma+i\omega)^k+s(n,k)(\gamma-i\omega)^k\right]\\
&=\sum_{k=0}^ns(n,k)r^k\cos(k\theta),
\end{align*}
where $(\cdot)_n$ is the falling factorial, $s(\cdot,\cdot)$ is the Stirling number of the first kind, and we use the substitution $\gamma=r\cos\theta$ and $\omega=r\sin\theta$. Alternatively, without the trig substitution, we can get
$$a_n=\sum_{k=0}^n s(n,k)\sum_{m=0}^k {k\choose m}\gamma^{k-m}\omega^m\cos\left(\frac{m\pi}{2}\right)$$
