The application of the finite element method to problems in neutron diffusion in space, energy, and time is studied. The use of piecewise polynomials with a variational form of the diffusion equation leads to algebraic systems of equations with characteristics similar to the usual finite difference equations. In Part I, a theoretical analysis of the finite element method, with Hermite polynomials, is presented and rigorous error bounds for the approximate solution are developed. In Part II, numerical studies are presented for problems in space and time. The results confirm the theoretical analysis and indicate the power of the method for diffusion problems.