In this paper, we consider the 3D problem of laser-induced semiconductor plasma generation under the action of the optical pulse, which is governed by the set of coupled time-dependent non-linear PDEs involving the Poisson equation with Neumann boundary conditions. The main feature of this problem is the non-linear feedback between the Poisson equation with respect to induced electric field potential and the reaction-diffusion-convection-type equation with respect to free electron concentration and accounting for electron mobility (convection’s term). Herein, we focus on the choice of the numerical method for the Poisson equation solution with inhomogeneous Neumann boundary conditions. Despite the ubiquitous application of such a direct method as the Fast Fourier Transform for solving an elliptic problem in simple spatial domains, we demonstrate that applying a direct method for solving the problem under consideration results in a solution distortion. The reason for the Neumann problem’s solvability condition violation is the computational error’s accumulation. In contrast, applying an iterative method allows us to provide finite-difference scheme conservativeness, asymptotic stability, and high computation accuracy. For the iteration technique, we apply both an implicit alternating direction method and a new three-stage iteration process. The presented computer simulation results confirm the advantages of using iterative methods.