The rotational velocities of stars in galaxies indicate that either Newtonian gravity breaks down on this scale, or that there is extra mass in the universe that has not been observed. The former is known as Modified Newtonian Dynamics (MOND), whereas the latter option is known as the dark matter paradigm. In this thesis a version of MOND is considered where the field equation for gravity is modified and becomes nonlinear. To simulate MOND, new N-body codes are needed due to the field equation not being linear anymore and in this thesis a new particle mesh method is developed. This new particle mesh method follows the same steps as the normal particle mesh method, except at the step where the acceleration field on the grid is calculated. This step is replaced by an iterative method to find the MONDian acceleration field on the grid. The code was then tested in the deep MOND regime using four cases of which analytical solutions of the MOND Poisson equation are known. These cases are: linear motion, the two-body problem, a ring system and an isothermal sphere. Using these test cases, numerous conclusions were drawn. From the linear motion test it was found that in the code bodies can interact with themselves, resulting in non-physical accelerations. The two-body test showed that the magnitude of the relative error in the force as calculated by the particle mesh code is roughly on the same order as the same error when calculating the Newtonian force using a particle mesh code. Furthermore, all tests showed that energy is conserved quite well, even though momentum need not be conserved due to the self interactions, and that generally the conservation of angular momentum is violated by the code. The tests also showed that the analytical formulae that were derived for these systems gave largely the same predictions as the algorithm, hence verifying both the formulae as the algorithm. The code could be sped up by writing it in C, the accuracy can be improved by reducing self-interactions and it can be extended by adding a method to handle close encounters. Some applications of the code are given, including: simulating the solar system, wide binary stars, tidal dwarf galaxies and galaxy clusters.