We consider a system of nonlinear wave equations with constraints that arises from the Einstein equations of general relativity and describes the geometry of the so-called Gowdy symmetric spacetimes on T3. We introduce two numerical methods, which are based on pseudo-spectral approximation. The first approach relies on marching in the future time-like direction and toward the coordinate singularity t=0. The second approach is designed from asymptotic formulas that are available near this singularity; it evolves the solutions in the past timelike direction from final data given at t=0. This backward method relies a novel nonlinear transformation, which allows us to reduce the nonlinear source terms to simple quadratic products of the unknown variables. Numerical experiments are presented in various regimes, including cases where spiky structures are observed as the coordinate singularity is approached. The proposed backward strategy leads to a robust numerical method which allows us to accurately simulate the long-time behavior of a large class of Gowdy spacetimes.