  We introduce a new class of singular partial differential equations, referred
to as the second-order hyperbolic Fuchsian systems, and we investigate the
associated initial value problem when data are imposed on the singularity.
First, we establish a general existence theory of solutions with asymptotic
behavior prescribed on the singularity, which relies on a new approximation
scheme, suitable also for numerical purposes. Second, this theory is applied to
the (vacuum) Einstein equations for Gowdy spacetimes, and allows us to recover,
by more direct arguments, well-posedness results established earlier by Rendall
and collaborators. Another main contribution in this paper is the proposed
approximation scheme, which we refer to as the Fuchsian numerical algorithm and
is shown to provide highly accurate numerical approximations to the singular
initial value problem. For the class of Gowdy spacetimes, the numerical
experiments presented here show the interest and efficiency of the proposed
method and demonstrate the existence of a class of Gowdy spacetimes containing
a smooth, incomplete, and non-compact Cauchy horizon.
