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.