Data from the Voyager probes and the Interstellar Boundary Explorer have revealed the importance of pick-up ions (PUIs) in understanding the character and behavior of the outer heliosphere, the region of interaction between the solar wind and the interstellar medium. In the outer heliosphere PUIs carry a large fraction of the thermal pressure, which effects the nature of the termination shock, and they are a dominate component of pressure in the heliosheath. This paper describes the development of a new multi-ion, multi-fluid 3-D magnetohydrodynamic model of the outer heliosphere. This model has the added capability of tracking the individual fluid properties of multiple ion populations. For this initial study two ion populations are modeled: the thermal solar wind ions and PUIs produced in the supersonic solar wind. The model also includes 4 neutral fluids that interact through charge-exchange with the ion fluids. The new multi-ion simulation reproduces the significant heating of PUIs at the termination shock, as inferred from Voyager observations, and provides properties of PUIs in the 3-D heliosheath. The thinning of the heliosheath due to the loss of thermal energy in the heliosheath from PUI and neutral interaction is also quantified. In future work the multi-ion, multi-fluid model will be used to simulate energetic neutral atom (ENA) maps for comparison with the Interstellar Boundary Explorer, particularly at PUI energies of less than 1 keV.