  The solution of systems of non-autonomous linear ordinary differential
equations is crucial in a variety of applications, such us nuclear magnetic
resonance spectroscopy. A new method with spectral accuracy has been recently
introduced in the scalar case. The method is based on a product that
generalizes the convolution. In this work, we show that it is possible to
extend the method to solve systems of non-autonomous linear ordinary
differential equations (ODEs). In this new approach, the ODE solution can be
expressed through a linear system that can be equivalently rewritten as a
matrix equation. Numerical examples illustrate the method's efficacy and the
low-rank property of the matrix equation solution.
