  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. As an application, the
method is used to compute rigorously stable and unstable bundles of periodic
orbits of vector fields. Examples are given in the context of the Lorenz
equations and the $\zeta^3$-model.
