  We introduce a new method with spectral accuracy to solve linear
non-autonomous ordinary differential equations (ODEs) of the kind $
\frac{d}{dt}\tilde{u}(t) = \tilde{f}(t) \tilde{u}(t)$, $\tilde{u}(-1)=1$, with
$\tilde{f}(t)$ an analytic function. The method is based on a new expression
for the solution $\tilde{u}(t)$ given in terms of a convolution-like operation,
the $\star$-product. This expression is represented in a finite Legendre
polynomial basis translating the initial problem into a matrix problem. An
efficient procedure is proposed to approximate the Legendre coefficients of
$\tilde{u}(t)$ and its truncation error is analyzed. We show the effectiveness
of the proposed procedure through some numerical experiments. The method can be
easily generalized to solve systems of linear ODEs.
