In computational relativity, critical behaviour near the black hole threshold has been studied numerically for several models in the last decade. In this paper we present a spatial Galerkin method, suitable for finding numerical solutions of the Einstein-Dirac equations in spherically symmetric spacetime (in polar/areal coordinates). The method features exact conservation of the total electric charge and allows for a spatial mesh adaption based on physical arclength. Numerical experiments confirm excellent robustness and convergence properties of our approach. Hence, this new algorithm is well suited for studying critical behaviour.