A new method for solving non-autonomous ordinary differential equations is
proposed, the method achieves spectral accuracy. It is based on a new result
which expresses the solution of such ODEs as an element in the so called
⋆-algebra. This algebra is equipped with a product, the ⋆-product,
which is the integral over the usual product of two bivariate distributions.
Expanding the bivariate distributions in bases of Legendre polynomials leads to
a discretization of the ⋆-product and this allows for the solution to be
approximated by a vector that is obtained by solving a linear system of
equations. The effectiveness of this approach is illustrated with numerical
experiments