I am interested in solving linear recurrences of the form $$a_{n+1}=\sum_{i=0}^K n^i X_i + \sum_{i=0}^L n^i Y_i a_n \tag{1}$$ where the $Y_i$ are $N\times N$ matrices, and the $X_i$ and $a_n$ are $N\times 1$ column vectors. (I am specifically interested in the $K=1,L=2,N=3$ case).
In the scalar ($N=1$) case, there is literature on how to find a closed-form expression for $a_n$ in terms of hypergeometric functions (Petkovšek's algorithm), and indeed computer algebra software returns an explicit solution for (1) in the scalar $N=1,K=1,L=2$ case.
Are there known techniques for solving (1) in the matrix case?