In this thesis we present a series of studies in numerical relativity investigating stability, hyperbolicity and critical phenomena. The first part of our work is dedicated to the study of d-stars, hypothetical objects consisting of a boson star and global monopole minimally or nonminimally coupled to the general relativistic gravitational field. The space of solutions for these systems is large and, for a wide range of coupling parameters, exhibits ground state solutions with asymptotic shells of bosonic matter. After demonstrating the existence of these stationary solutions, we turn our attention to their stability through a combination of perturbation theory and dynamical simulation. In doing so, we demonstrate that the novel solutions we have found, as well as the highly compact solutions investigated by previous authors, appear to be generically unstable to radial perturbations. As such, we find that d-stars are poor candidates for astrophysically relevant black hole mimickers. Generalizing from the stability of solutions to the stability of methods, we introduce a novel formulation of numerical relativity which we refer to as reference metric covariant and conformal Z4 (RCCZ4). Like its Z4 namesake, RCCZ4 promotes the 3+1 Hamiltonian and momentum constraints to dynamical degrees of freedom. Unlike Z4 however, RCCZ4 accomplishes this by coupling the constraints to a reference metric completely independent of the physical metric. Although we have only investigated RCCZ4 in the case of time independent Lorenzian reference metrics, the method may generalize to user specifiable reference metrics which could potentially confer additional beneficial properties. Even in this simple case, however, the performance of RCCZ4 is comparable to leading hyperbolic formulations. The final part of our thesis works towards developing superior understanding of strong field gravity through the investigation of gravitational collapse. We consider the system consisting of the electromagnetic and general relativistic gravitational fields and investigate the threshold of black hole formation in axisymmetry. Previous studies of this system have reported family dependent scaling phenomena as criticality is approached. Although we find good agreement with previous investigations of dipole-type initial data, our investigations of quadrupole-type initial data point towards universal scaling as opposed to family dependent scaling.