In this thesis, we investigate statistical properties of the redshift zero matter distribution in disformal gravity through N-body simulations. The disformal model studied here is a conformally coupled field in a symmetron potential, but with an additional exponential disformal term. First, important concepts about general relativity, modified gravity, and simulations are presented. The equation of motion for the disformal scalar field is found, along with the modified geodesics of dark matter particles in this model. We implement the resulting differential equations into the N-body code Isis/Ramses. The results of tests, which compare results from the disformal code to simulations done with earlier symmetron-only codes, are shown; and the code reproduces the symmetron results fairly well. Finally, we perform cosmological simulations with 256^3 particles, for five different parameter sets, for which the matter power spectrum and the halo mass function are shown. We also present the average field profile, the field oscillation amplitude, and the magnitude of the fifth forces around a massive halo. The conformally coupled symmetron increases both the power spectrum and the mass function noticeably, compared to LCDM. The main result of this study is that adding a strong disformal term can mask some of the increase in clustering, thereby bringing both the power spectrum and the mass function closer to general relativity. Furthermore, we found that the disformal term gives rise to oscillations of the scalar field in high density regions. Such oscillations can lead to increased magnitudes of the fifth forces, compared to the symmetron alone.