A methodology is presented for the systematic modal reduction of structural systems which contain quadratic and cubic nonlinearities in displacement. The procedure is based on the center manifold approach for describing individual nonlinear modes, but it has been extended to account for simultaneous motion within several chosen modal coordinates. Motions of the reduced system are constrained to lie on high-dimensional manifolds within the phase space of the original system. Polynomial approximations of these manifolds are obtained through third order for arbitrary system parameters. Algorithms have been developed for automation of this procedure, and they are applied to an example system. Free and forced responses of the reduced system are discussed and compared to responses reduced through simple modal truncation. A more rigorous treatment of harmonic forcing is proposed, which will allow for the production of high-dimensional, time-dependent manifolds through a simple adaptation of the unforced procedure.