The rotational velocities of stars at the edges of galaxies are higher than predicted by Newtonian dynamics. One theory that is capable of explaining this anomaly is called MOdified Newtonian Dynamics
(MOND). To further test this MOND theory, the wide binary test has been suggested. In the wide binary test, the stability of wide binary systems is tested. The MOND theory deviates significantly from the
Newtonian theory when accelerations are smaller than a0 = 1.2·10−10 m/s2 , which is present in wide binary systems. The MOND theory is non-linear, which introduces a non-trivial effect of the external field
(EF) of the Milky Way on the wide binary orbit. To perform the wide binary test, an iterative particle mesh N-body code was used with initial conditions taken from the Gaia mission database. Additionally, the EF was modeled in three ways, where the Newton potential adjustment model was found to be the most realistic. The results include nine simulations, split up into three models and three different regimes. The models include the Newton model, the MOND without EF model, and the MOND with EF model. The regimes consist of the Newtonian regime (ainternal > a0), the transition regime (ainternal ≈ a0), and the deep MOND regime (ainternal < a0). In the Newton model, only the Newtonian regime gave stable orbits. The MOND without EF model gave stable orbits for all regimes, and the MOND with EF model gave stable orbits only for the Newtonian regime. The MOND without EF gave stable wide binary orbits, whereas the other models did not. If additional measurements detect these wide binary systems, then the MOND without EF model gains credibility. The orbits of the Newton model and the MOND with EF model were similar, with only a small difference. This is because the strength of the external field is 1.6 · a0, which makes the external field dominant over the internal accelerations, making the simulation Newtonian-like. To determine which model is correct, additional long-term measurements of Gaia are required. Further research recommendations include improving the performance and accuracy of the code. This is done by better implementation of the fast Fourier transform, translating the code into a faster programming language than Python, using adaptive time steps, and reducing self-interaction. The most important extension to the model is to introduce the galactic tidal effect.