The mathematical model of investigated problem [1, 2] represents the system of nonlinear differential equations with respect to unknown functions that are electrostatic potential, electron concentration and holes concentration. The problem is significantly complicated by the fact that the boundary conditions are inhomogeneous: on a part of the boundary the Dirichlet conditions are specified, and on the other part – the Neumann conditions. One of the most important points in the numerical solution of this problem is the problem of choosing a suitable method for linearizing the given differential equations. For this purpose, we use two different approaches: Newton’s and Gummel’s algorithms [2]. Then the obtained linear system is solved numerically by applying the conjugate and bi-conjugate gradient methods. Numerical experiments were carried out on the basis of these two techniques. The obtained numerical solutions made it possible to compare the considered methods and to draw preliminary conclusions about the possibilities of their application.