The study of binary neutron star coalescence, one the most energetic classes of events in the universe, requires calculating the complicated interactions of strong gravity, relativistic fluids, and magnetic fields. The Spectral Einstein Code provides a framework for simulating the inspiral and merger of black holes and neutron stars, but its ability to model the behavior of binary neutron stars and magnetic fields is a recent development. This work describes the implementation of an initial data solver for neutron star binaries, a magnetohydrodynamics module for neutron star and accretion disk evolutions, a pair of basis functions wellsuited to spectral representations of neutron star spacetimes, and a selection of other improvements to this research code. It also presents the results of early investigations using these new capabilities, including the effects of magnetic fields on shear instabilities in differentially rotating neutron stars. Such stars may be formed from core-collapse supernovae or low-mass binary neutron star mergers, and fluid instabilities in galactic sources can produce gravitational waves observable by detectors in the near future. We find that strong magnetic fields are capable of suppressing a shear instability, but they also trigger magnetic instabilities whose effects may be just as observable as the original signal.