Plasma fluid modeling equations such as continuity equations for charged species with drift-diffusion approximation electron energy equation, and Poisson's equation are solved by using a new paradigm. Resulting discretized equations are solved jointly by the Newton-Krylov-Schwarz (NKS) scheme by means of a parallelized toolkit called PETSc. All model equations are nondimensionalized and are discretized using fully implicit finite-difference method with the Scharfetter-Gummel for the fluxes. Thermal flux is considered for electrons, while both thermal and drift fluxes are considered for ions at the electrodes. Both LU and ILU linear equation solvers are tested in solving the preconditioned matrix equation in each subdomain. Results show that with the present test case the combination of Block Jacobi preconditioner with ILU and BiCG-Stab performs the best with 92.8% of parallel efficiency at 28 processors.