A hydrodynamic plasma simulation of the ionized component of the global heliosphere has been coupled to a nonstationary Boltzmann particle-mesh method which describes the neutral hydrogen present in the same region. This coupling is a numerical model of the physical charge-exchange process through which the two particle species interact. Numerical self-consistency is achieved by iterating the distributions of the two particle species between the two numerical schemes until a steady state is achieved. This approach will contribute to the understanding of the solar wind—LISM interaction. ©1999 American Institute of Physics.