  This paper develops a computational method for studying stable/unstable
manifolds attached to periodic orbits of differential equations. The method
uses high order Chebyshev-Taylor series approximations in conjunction with the
parameterization method -- a general functional analytic framework for
invariant manifolds. The parameterization method can follow folds in the
embedding, recovers the dynamics on the manifold through a simple conjugacy,
and admits a natural notion of a-posteriori error analysis. The key to the
approach is the derivation of a recursive system of linear differential
equations describing the Taylor coefficients of the invariant manifold. We find
periodic solutions of these equations by solving a coupled collection of
boundary value problems with Chebyshev spectral methods. We discuss the
performance of the method for the Lorenz system, and for circular restricted
three and four body problems. We also illustrate the use of the method as a
tool for computing cycle-to-cycle connecting orbits.
