We computationally explore the dynamics of superconductivity near the superheating field in two ways. First, we use a finite element method to solve the time-dependent Ginzburg-Landau equations of superconductivity. We present a novel way to evaluate the superheating field Hsh and the critical mode that leads to vortex nucleation using saddle-node bifurcation theory. We simulate how surface roughness, grain boundaries, and islands of deficient Sn change those results in 2 and 3 spatial dimensions. We study how AC magnetic fields and heat waves impact vortex movement. Second, we use automatic differentiation to abstract away the details of deriving the equations of motion and stability for Ginzburg-Landau and Eilenberger theory. We present calculations of Hsh and the critical wavenumber using linear stability analysis.