The interaction between plasma, neutrals and surfaces in fusion reactors is of immense importance for heat and particle control, especially for the next generation of devices. Heat loads of are expected for steady state operation at ITER and up to in slow transient situations. To study the complex physics appearing between the plasma and the divertor material, as well as techniques for heat flux mitigation, plasma linear devices are employed. Magnum-PSI, located at DIFFER, can reproduce the heat and particle loads expected at ITER. However, due to the complexity of the plasma–wall interaction, numerical models are required to better understand the experiments and to extrapolate the results to a tokamak divertor configuration. For tokamak geometries, SOLPS-ITER (formerly known as B2.5-Eirene) is employed to solve the plasma and neutral distribution in a coupled way. However, the utilization of this code for linear devices is not straightforward. Thus, a neutral module was developed with linear devices in mind, named Eunomia. Nevertheless, there is still a relevant interest in using SOLPS-ITER with linear devices, as it allows to easily transfer knowledge about relevant atomic and molecular processes close to the surface and the effect of different mitigation techniques. This work presents a systematic comparison between the two neutral modules, Eirene and Eunomia, in stand-alone and coupled runs. Special attention is paid to the implementation of plasma-neutral interactions, in which both codes diverge significantly. The sources of particles and energy that are used by B2.5 in a coupled run are analysed. Significant differences in the implementation of electron impact ionization, molecular assisted recombination and proton–molecule elastic collisions lead to disparate sources of particles and energy and, in some cases, differences in the distribution of neutrals achieved by each code. Moreover, a double counting in proton–atom collisions was identified in Eunomia as a result of this analysis, artificially increasing the plasma-neutral sink of energy. This would lead to different plasma evolutions in coupled runs. Nevertheless, additional-free parameters in both coupled code suites leave sufficient freedom to match experimental data. Additional data would be required to further constrain these parameters and the coupled solution.