One of the primary concerns in galaxy evolution is how galaxies form their stars: what keeps that star formation going over cosmic time, and what causes it to stop in a processes called "quenching". Galaxies with mass similar to our own Milky Way occupy a sweet spot between abundance and brightness that makes them easy to find in the sky, and such galaxies also populate a transitionary regime in behavior that make them interesting for studying galaxy evolution. Numerical modeling-from semi-analytic models to numerical simulations-are valuable tools for understanding the multiple intersecting physical processes that drive galaxy evolution. These processes act both within and around individual galaxies such that numerical models must necessarily encompass a range of spatial and temporal scales. Multiple approaches are commonly used in order for this modeling to be physically insightful. In this dissertation I will present my efforts to unravel the mechanisms of galaxy evolution affect Milky Way-like galaxies using a variety of numerical models.Addressing the issue of what causes galaxies to stop forming stars, I first investigate an unusual population of galaxies called the "break BRDs" (Tuttle and Tonnesen 2020). Within the dominant framework for galaxy quenching, galaxies first stop forming stars in their centers and later in their outskirts. This is the "inside-out" quenching paradigm. The break BRD galaxies possess observational markers that run counter to this narrative. We used the Illustris TNG cosmological simulation(Pillepich et al. 2018b) to find a set of simulated galaxies that are analogous to the observed breakBRDs in order to better understand their evolution. We found that the breakBRD analogues are galaxies that ultimately become fully quenched, but found no clear cause for the "outside-in" modality. This is not the dominant channel for quenching in the IllustrisTNG simulation, but roughly 10% of quiescent galaxies with 10