We describe an implementation of density functional theory that is formulated fully in configuration space, where all wave functions, densities, and potentials are represented on a grid. Central to the method is a fourth-order factorization of the evolution operator for the Kohn-Sham Hamiltonian. Special attention is paid to nonlocal pseudopotentials of the Kleinman-Bylander type, which are necessary for a quantitatively accurate description of molecules, clusters, and solids. It is shown that the fourth-order factorization improves the computational efficiency of the method by about an order of magnitude compared with second-order schemes. We use the Ono-Hirose filtering method to reduce the resolution of the grid used for representing the wave functions. Some care is needed to maintain the fourth-order convergence using filtered projectors, and the necessary precautions are discussed. We apply the method to an isolated carbon atom as well as the carbon-monoxide molecule, the benzene molecule, and the buckminsterfullerene cluster, obtaining quantitative agreement with previous results. The convergence of the method with respect to time step, grid resolution, and filtering method is discussed in detail.