In the high energy limit of hadron collisions, the evolution of the gluon density in the longitudinal momentum fraction can be deduced from the Balitsky hierarchy of equations or, equivalently, from the nonlinear Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation. The solutions of the latter can be studied numerically by using its reformulation in terms of a Langevin equation. In this paper, we present a comprehensive study of systematic effects associated with the numerical framework, in particular the ones related to the inclusion of the running coupling. We consider three proposed ways in which the running of the coupling constant can be included: "square root" and "noise" prescriptions and the recent proposal by Hatta and Iancu. We implement them both in position and momentum spaces and we investigate and quantify the differences in the resulting evolved gluon distributions. We find that the systematic differences associated with the implementation technicalities can be of a similar magnitude as differences in running coupling prescriptions in some cases, or much smaller in other cases.