Computing the stable, unstable and center manifolds of a given dynamical system gives one much insight into the properties and evolution of the system. For this reason, several techniques for their determination (both analytical and numerical) have been proposed over the years. In the particular case of center manifolds, numerical techniques based on Taylor polynomials provide reasonable approximations, but only in a neighborhood of an equilibrium point or a periodic orbit. The aim of the present paper consists in approximating the center manifold by certain functions obeying a system of partial differential equations. These functions are subsequently determined by applying Chebyshev spectral methods. As a result, the validity of the approximations is extended to larger regions in the phase space. The technique is described in detail and illustrated for several dynamical systems, both Hamiltonian and non-Hamiltonian, with one-, two- and four-dimensional center manifolds. Although the technique shows accurate results for the examples considered, computing higher-dimensional center manifolds still presents challenging computational issues.