The electron component of an ultracold neutral plasma (UCP) is modeled based on a scalable method using a self-consistently determined mean-field approximation. Representative sampling of discrete electrons within the UCP are used to project the electron spatial distribution onto an expansion of orthogonal basis functions. A collision operator acting on the sample electrons is employed in order to drive the distribution toward thermal equilibrium. These equilibrium distributions can be determined for non-zero electron temperatures even in the presence of spherical symmetry-breaking applied electric fields. This is useful for predicting key macroscopic UCP parameters, such as the depth of the electrons' confining potential. Dynamics such as electron oscillations in UCPs with non-uniform density distributions can also be treated by this model.