Recurrence relation with non-constant coefficients and modular arithmetic

Question

Let $a_n$ be a linear homogenous recurrence relation with constant coefficients. One can use the matrix form of the recurrence and then it is easy to determine $a_k$ mod $p$ for all $k$ and any prime $p$. This method would work independent of the order of the recurrence relation.
Now, assume that we have another linear homogenous recurrence relation but with non-constant coefficients, say $b_n$. Is there a method to determine $b_k$ mod $p$ with $k$ and $p$ being the same as above? Here, a general answer (independent of the order of the recurrence) would be appreciated. In case a general answer isn't possible, you can assume that $b_n$ is a $3^{rd}$ order recurrence relation of the form $f(n)b_n = g(n)b_{n-1} + h(n)b_{n-2} + i(n)b_{n-3}$ for distinct real valued functions $f$, $g$, $h$ and $i$.
If there does exist a method, would it be computationally efficient?