The hydromagnetic dynamo Z-model represents a nonlinear dynamic system. Its steady-state solution is derived by step-by-step integration of parabolic partial differential equations of the second order with the use of the finite difference method. Until now, two methods have been used: the semi-implicit method in which the θ-diffusion was carried out implicitly while the r-diffusion explicitly, and the implicit method in which the complete diffusion was carried out implicitly. In the present contribution, a combined semi-implicit method is suggested which reflects not only the singularity at the coordinate system origin but also the decreasing mesh size near the core-mantle boundary. This procedure preserves the advantage of semi-implicit methods and, simultaneously, increases the stability in the most critical boundary layer.