In this paper, a new rigorous numerical method to compute fundamental matrix
solutions of non-autonomous linear differential equations with periodic
coefficients is introduced. Decomposing the fundamental matrix solutions
$\Phi(t)$ by their Floquet normal forms, that is as product of real periodic
and exponential matrices $\Phi(t)=Q(t)e^{Rt}$, one solves simultaneously for
$R$ and for the Fourier coefficients of $Q$ via a fixed point argument in a
suitable Banach space of rapidly decaying coefficients.