A handful of Milgromian \(N\) -body codes were developed before por to handle MOND computations, and these have been applied to various problems. The first multi-grid, Milgromian \(N\) -body code was developed by [1]} to investigate the stability of disc galaxies. This was later extended to simulate how they might warp due to the external field effect [2]}. Another \(N\) -body solver which implemented the AQUAL formulation of MOND was developed and used to study the evolution of spiral galaxies using pure stellar discs [3]}. Gas dynamics was later included using a sticky particle scheme at low resolution [4]}. n-mody was developed to solve the Milgromian Poisson equation in spherical coordinates [5]} and used to investigate dynamical friction [6]}, orbit instabilities [7]}, and stellar kinematics [8]}, [9]}, [10]}. Milgromian \(N\) -body codes tailored to cosmological simulations have also been developed [11]}, [12]}, [13]}, [14]}. Another \(N\) -body solver called raymond [15]} was developed to implement both the AQUAL [16]} and QUMOND [17]} formulations of MOND. raymond has been applied to cosmological [18]}, galaxy cluster [19]}, and other problems. However, not all the aforementioned \(N\) -body codes can be applied to generic scenarios simultaneously involving particles, gas dynamics, and star formation.
