As the number of solar system small body multiple systems (including binaries) grows in number and in observational baseline, advanced dynamical modelling will reveal physical, orbital, and spin properties that are currently poorly understood. For example, among Kuiper Belt Objects (KBOs) beyond Neptune, only a few have well determined obliquities, but observations of orbital precession due to non-spherical shapes can identify spin pole directions and obliquities. The ~50 KBO Binaries with known orbital solutions have only been studied assuming Keplerian motion, but evidence suggests that several of these are beginning to show observationally significant non-Keplerian effects. To study these effects, we have developed and validated an n-quadrupole integrator SPINNY (SPIN + N-bodY) which can include all the dominant non-Keplerian effects relevant within observational precision. SPINNY has been combined with advanced Bayesian model fitting for KBO astrometry in a new modeling tool called MultiMoon. We present MultiMoon along with some potential use cases to demonstrate its value. We also present results on an analysis of the 47171 Lempo (1996 TC36) KBO trinary system using MultiMoon.