For the design of organic semiconductor devices such as organic light-emitting devices and solar cells, it is of crucial importance to solve the underlying charge transport equations efficiently and accurately. Only a fast and robust solver allows the use of fitting algorithms for parameter extraction and variation. Introducing appropriate models for organic semiconductors that account for the disordered nature of hopping transport leads to increasingly nonlinear and more strongly coupled equations. The solution procedures we present in this study offer a versatile, robust, and efficient means of simulating organic semiconductor devices. They allow for the direct solution of the steady-state drift-diffusion problem. We demonstrate that the numerical methods perform well in combination with advanced physical transport models such as energetic Gaussian disorder, density-dependent and field-dependent mobilities, the generalized Einstein diffusion, traps, and its consistent charge injection model.