Given a three-term recurrence relation such as
\begin{equation} a_nX^n+b_nX^{n+1}+c_nX^{n-1}=\zeta_n \delta_{n,0} \end{equation} it can always be written in continued-fraction form, and if you are lucky, it can be resumed (using Gauss's continued fraction). But this only works if these coefficients are just numbers.
Unfortunately, what I have is
\begin{equation} a_nZ^n+b_nZ^{n+1}+c_nZ^{n-1}=\zeta_n \delta_{n,0} \end{equation} where $Z^n=(X^n,Y^n)$ is vector and $a_n$, $b_n$, and $c_n$ are $2\times 2$ matrices. Therefore, I cannot use Gauss's continued fraction.
By working it out, I can express it as two-coupled continued-fraction:
\begin{equation} \begin{split} &a^{1}_nX^n+b^{1}_nX^{n+1}+c^{1}_nX^{n-1}=d^{1}_nY^n+e^{1}_nY^{n+1}+f^{1}_nY^{n-1}\\ &a^{2}_nX^n+b^{2}_nX^{n+1}+c^{2}_nX^{n-1}=d^{2}_nY^n+e^{2}_nY^{n+1}+f^{2}_nY^{n-1}\\ \end{split} \end{equation} where now the coefficients are just numbers and they are not the same, $\eta^{1}_n\neq\eta_n^{2}$ (I do not write them down as they are complicated expressions).
Is there any procedure I could follow to isolate the variables and get a three-term recursive relation for each $X$ and $Y$ independently? Otherwise, what would be the procedure to solve this coupled system?
Example: Suppose the following toy model \begin{equation} \begin{split} &X^n+2X^{n+1}+5X^{n-1}=2Y^n+4 Y^{n+1}+10Y^{n-1}+2\\ &X^n+3X^{n+1}+11X^{n-1}=5Y^n+8 Y^{n+1}+7Y^{n-1}+4 \end{split} \end{equation}
Solution for the toy model: The main procedure is to use the generating functions $X(z)=\sum x_n z^n$ and $Y(z)=\sum y_n z^n$. By plugging those expressions, one gets to coupled equations as
\begin{equation} \begin{split} &2\frac{X(z)-x_0-z x_1}{z^2}+\frac{X(z)-x0}{z}+5X(z)=4\frac{Y(z)-y_0-z y_1}{z^2}+2\frac{Y(z)-y0}{z}+10X(z)+2\frac{1}{1-z}\\ &3\frac{X(z)-x_0-z x_1}{z^2}+\frac{X(z)-x0}{z}+11X(z)=8\frac{Y(z)-y_0-z y_1}{z^2}+5\frac{Y(z)-y0}{z}+7X(z)+4\frac{1}{1-z} \end{split} \end{equation}
Those equations can now be solved for $X(z)$ and $Y(z)$ in terms of partial fractions and those can be used to find the isolated recurrence relations.
While this works for my toy model where the coefficients are constants, it is still unclear if this would work for more general, and $n-$dependent coefficients.
