The 2D scrape-off-layer turbulence code (SOLT) is extended to include neutral-plasma interactions. A Boltzmann equation is derived for the evolution of the bi-normally-averaged neutral distribution function, G(x,v x ,t), in the radial dimension, and this evolution is included in the new code (nSOLT). Neutral-plasma interactions are mediated by charge-exchange (CX) and ionization rates based on poloidally-averaged plasma density and temperature. Good agreement is obtained between asymptotically stationary neutral density profiles from nSOLT simulations and those previously obtained from the Monte Carlo neutral transport code DEGAS 2, for timeaveraged NSTX H-mode plasma profiles. The sensitivity of the nSOLT neutral profiles to atomic physics parameters, with and without CX physics, is included in the comparison. In addition, nSOLT simulations that evolve the plasma in 1D, using radial diffusion as a proxy for turbulent (blob) transport, illustrate the convergence to a self-consistent neutral-plasma equilibrium sustained by a neutral source at the far-SOL boundary and plasma heating in the core; equilibria consistent with typical NSTX Ohmic L-mode plasmas are described.