La e.g.f.'s$\sum_{n=0}^\infty r_n t^n/n!$, para las dos soluciones linealmente independientes
son
$ \exp(b t^2/2 + a t)$ $\exp(b t^2/2 + a t)\; \text{erf}(\sqrt{b/2} t + a/\sqrt{2b})$.
De la primera, obtenemos
$$ r_n = n! \sum_{k=0}^{\lfloor n/2 \rfloor} \dfrac{ (b/2)^k a^{n-2k}}{k! (n-2k)!}$$
EDIT: he Aquí una pequeña explicación. La e.g.f. (exponencial de generación de función) de una secuencia $r_n$ es la función de $g(t) = \sum_{n=0}^\infty r_n t^n/n!$.
Esto tiene la propiedad $g'(t) = \sum_{n=0}^\infty r_{n+1} t^n/n!$,
$g''(t) = \sum_{n=0}^\infty r_{n+2} t^n/n!$, etc., mientras
$$ \sum_{n=0}^\infty n\; r_n \dfrac{t^n}{n!} = \sum_{m=0}^\infty r_{m+1}\dfrac{ t^{m+1}}{m!} = t g'(t)$$
Ahora escriba su recurrencia como
$$ r_{n+2} = a \; r_{n+1} + b\; (n+1) r_{n} = a \; r_{n+1} + b \; r_n + b\; n r_n$$
Multiplicar cada término por $t^n/n!$ y suma. Tenemos
$$ g''(t) = a g'(t) + b g(t) + b t g'(t) $$
y dos soluciones linealmente independientes para esta ecuación diferencial se
$$ g(t) = \exp(b t^2/2 + a t) \ \text{and} \ g(t) = \exp(b t^2/2 + a t)\; \text{erf}\left(\sqrt{b/2}\; t + a/\sqrt{2b}\right)$$
Si $a$ $b$ son matrices que no conmutan, las cosas no son tan sencillas:
ya no es cierto que $\dfrac{d}{dt} \exp(b t^2/2 + a t) = ( bt + a) \exp(b t^2/2 + a t)$. Incluso para el $2 \times 2$ caso no creo que te voy a obtener de forma cerrada soluciones.