Fluid descriptions of plasmas, which are usually applied to a collisional plasma, can only be justified for very small Coulomb Knudsen numbers. However, the scrape-off layer (SOL) plasmas of experimental magnetic confinement fusion devices tend to have operational regimes characterized by a Coulomb Knudsen number around 0.1. In interesting detached regimes of an SOL plasma in a tokamak, when the plasma detaches from the limiters or divertors, this number may increase along with the local plasma gradients. Plasma gradients are also known to increase (and thus drive non-local effects) in inertial confinement fusion. Neutrals, which are being produced owing to plasma recombination at the plasma–divertor interface, may be in a mixed collisional regime as well. Thus simultaneous kinetic treatments of plasma and neutral particles with self-consistent evaluation of boundary conditions at the material walls are required. We present a physical model and a numerical scheme, and discuss results of purely kinetic simulations of plasmas and neutrals for actual conditions in the Alcator C-Mod and Tokamak-de-Varennes experimental tokamaks. Results for both steady-state and transient regimes of SOL plasma flow are presented. Our approach, unlike particle-in-cell and Monte Carlo methods, is free from statistical noise.