Simulations of the evolution of the universe under gravity are essential to our understanding of modern cosmology. These $N$-body simulations can contain trillions of particles and are run on some of the most powerful supercomputers in the world. $N$-body simulations are necessary to understand the data from large galaxy surveys that are mapping the universe in higher detail than ever before. This thesis explores the accuracy of $N$-body simulation methods and their applications towards the Dark Energy Spectroscopic Instrument (DESI) Survey. Firstly, the accuracy of the gravity scheme in the Smoothed Particle Hydrodynamics With Inter-dependent Fine-grained Tasking (SWIFT) simulation code is tested. These tests inform us on the limitations of running SWIFT with high values of the opening angle parameter $theta$. Additionally, an error in the large-scale clustering in SWIFT simulations was found and fixed. Next, a comparison between several simulation codes is presented. In this comparison, simulations were run from identical initial conditions and the level of discrepancy caused by choice of simulation code is measured. The systematic errors caused by choice of code are compared to the statistical errors in a DESI Survey volume. We find that the matter power spectra from independent codes agree to within 1% for $k