Generating functions seem to be a powerful tool in discrete mathematics for solving differential equations and recurrence relations. I've been trying to figure out if these methods can be expanded to differential equations that involve matrices, such as Schrodinger's equation. Is there anything that prevents a solution to these types of differential equations being written using generating functions? For example, given a solution to Schrodinger's equation as
$$ U = Te^{-i \hbar \int_{0}^{t'} H(t) \,dt} $$
Where H(t) is some time-dependent Hamiltonian matrix and T is the time-ordering operator. This form seems very reminiscent to exponential generating functions. For example, the generating function for the Bessel functions is given by $$ e^{\frac{x}{2}(t-\frac{1}{t})}=\sum_{n=-\infty }^{\infty} J_n(x)t^n $$
So, if the time evolution of the Hamiltonian gave something similar to the generating function above, where x is now the matrix Hamiltonian, can the exponential be re-written in a form using the Bessel functions of matrix argument? I'd assume that the matrix argument of a Bessel function would be similar to an analytic function of matrix argument, but everything I find seems to write Bessel/Hypergeometric functions in terms of zonal polynomials, and the wikipedia page for hypergeometric functions of matrix argument even mention that these types of functions aren't similar to writing other functions in terms of matrix arguments. That doesn't make sense to me though since these special functions have these generating function relations. Any help would be appreciated if someone could point me towards the literature too.