We describe a rapidly converging algorithm for solving the Schrodinger equation with local potentials in real space. The algorithm is based on solving the Schrodinger equation in imaginary time by factorizing the evolution operator e(-epsilonH) to fourth order with purely positive coefficients. The wave functions \ psi (j)> and the associated energies extracted from the normalization factor e(j)(-epsilonE) converge as O(epsilon (4)). The energies computed directly from the expectation value, < psi (j)\H \ psi (j)>, converge as O(epsilon (8)). When compared to the existing second-order split operator method, our algorithm is at least a factor of 100 more efficient. We examine and compare four distinct fourth-order factorizations for solving the sech(2)(ax) potential in one dimension and conclude that all four algorithms converge well at large time steps, but one is more efficient. We also solve the Schrodinger equation in three dimensions for the lowest four eigenstates of the spherical analog of the same potential. We conclude that the algorithm is equally efficient in solving for the low-lying bound-state spectrum in three dimensions. In the case of a spherical jellium cluster with 20 electrons, our fourth-order algorithm allows the use of very large time steps, thus greatly speeding up the rate of convergence. This rapid convergence makes our scheme particularly useful for solving the Kohn-Sham equation of density-functional theory and the Gross-Pitaevskii equation for dilute Bose-Einstein condensates in arbitrary geometries.