A new method for investigating the interaction of the solar wind with the partially ionized local interstellar medium (LISM) is presented. The solar wind and the interstellar plasma are modeled using a two-dimensional (2-D) hydrodynamic numerical code. The plasma is coupled to the neutral hydrogen (of both interstellar and solar wind origin) via resonant charge exchange. To model the neutral H distribution, we use a nonstationary 2.5-D particle mesh method to solve the Boltzmann equation, which is coupled self-consistently to the interstellar and solar wind plasma. Numerical self-consistency is achieved by iterating the plasma and neutral H distributions between the two numerical schemes until a steady state is achieved. Results from three test applications are presented and discussed, including the first one-shock kinetic simulation. The simulations are able to reproduce the main features of the heliosphere such as shock structure, hydrogen wall, and heating, deceleration and filtration of neutral hydrogen. In addition, they enable the study and interpretation of the non-Maxwellian hydrogen distribution function. Traces of fast neutrals originating inside the termination shock and the heliosheath/heliotail region can be found far upstream of the outer heliosphere. The influence of different interstellar plasma boundary values on the heliosphere is highlighted in the comparison of two supersonic simulations and one subsonic simulation. In particular, by comparing the simulated energetic neutral atom (ENA) fluxes at 1 AU of the supersonic and subsonic models, it is found that the subsonic flux is significantly underabundant in the energy range 10--60 eV compared to the supersonic case. This may offer an important diagnostic for determining whether the heliosphere possesses a bow shock or not. ¿ 2000 American Geophysical Union