In this thesis, we describe numerical spherical collapse solutions to a ``modified gravity'' theory, Einstein dilaton Gauss-Bonnet (EdGB) gravity. Of the class of all known modified gravity theories, EdGB gravity has attracted recent attention due to speculations that the theory may have a classically well-posed initial value formulation and yet also exhibit stable scalarized black hole solutions (what makes this surprising is the plethora of black hole ``no hair theorems'', the assumptions behind which EdGB gravity manages to avoid). If EdGB gravity indeed possess these properties, it would be an ideal theory to perform model-selection tests against general relativity (GR) in binary black hole merger using gravitational waves. Furthermore, the theory is an important member of the so-called ``Horndeski theories'', which have been invoked to construct, e.g. nonsingular black hole and cosmological solutions, and to address the classical flatness and horizon problems of the early universe.

In constructing numerical solutions to EdGB gravity (without any approximations beyond the restriction to spherically symmetric configurations), we are able to carefully examine various claims made in the literature about EdGB gravity, perhaps most importantly whether or not the theory admits

a classically well-posed initial value problem. One conclusion of these studies has been, at least in spherical collapse, EdGB gravity can *dynamically* lose hyperbolicity, which shows EDGB gravity is fundamentally of ``mixed type''. Mixed type problems appear in earlier problems in mathematical physics, perhaps most notably in the problem of steady state, inviscid, compressible fluid flow. The loss of hyperbolicity and subsequent formation of ``elliptic regions'' outside of any sort of horizon implies that the theory violates cosmic censorship, broadly defined. Arguably it is clear that this result is gauge invariant, although we do not formulate a rigorous proof that this is so. In addition to discussing the hyperbolicity of EdGB gravity, we discuss several other interesting features to the numerical solutions, including the formation of scalarized black hole solutions in the theory, at least for certain parameter ranges for the theory and certain open sets of initial data.