An efficient way of evolving a solution to an ordinary differential equation is presented. A finite element method is used where we expand in a convenient local basis set of functions that enforce both function and first derivative continuity across the boundaries of each element. We also implement an adaptive step-size choice for each element that is based on a Taylor series expansion. This algorithm is used to solve for the eigenpairs corresponding to the one-dimensional soft Coulomb potential, 1/sqrt[x(2)+β(2)], which becomes numerically intractable (because of extreme stiffness) as the softening parameter (β) approaches zero. We are able to maintain near machine accuracy for β as low as β = 10(-8) using 16-digit precision calculations. Our numerical results provide insight into the controversial one-dimensional hydrogen atom, which is a limiting case of the soft Coulomb problem as β → 0.