Orbital systems of two or more Trans-Neptunian Objects (TNOs) are valuable for investigating formation processes of the early Solar System. Because TNOs exist below the diffraction limit of current telescopes, they are rendered as point spread functions (PSFs) in observational images. To correctly understand the orbital dynamics of TNO systems, the PSFs must be modeled according to telescope, camera, and observational parameters. These models can then be fit to images to produce the precise relative astrometry that is needed to fit the orbits. The software package nPSF was designed to fit modeled PSFs to images containing an n number of PSFs using Bayesian parameter inference methods, specifically the Markov Chain Monte Carlo process. The relative astrometry of the system can then be derived from the posterior distribution. We test the capabilities of nPSF by fitting images of 2PSF Trans-Neptunian Binaries and comparing our results to the published astrometry for those systems. We show that nPSF produces good astrometry in the tested systems. Additionally, we apply nPSF to the search of a tertiary object in various potential hierarchical triple systems. The results produced by nPSF do not indicate the detection of a third component in these systems.