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.