As semiconductor technology continues to evolve, numerical modeling of semiconductor devices becomes an indispensible tool for the prediction of device characteristics. The simple drift-diffusion model is still widely used, especially in the study of subthreshold behavior in MOSFETs. The numerical solution of these two equations offers difficulties in small devices and special methods are required for the case when dealing with 3D problems that demand large CPU times. In this work we investigate the convergence properties of the Bi-CGSTAB method. We find that this method shows superior convergence properties when compared to more commonly used ILU and SIP methods.