We discretize a unipolar electrothermal drift-diffusion model for organic semiconductor devices with Gauss–Fermi statistics and charge carrier mobilities having positive temperature feedback. We apply temperature dependent Ohmic contact boundary conditions for the electrostatic potential and use a finite volume based generalized Scharfetter-Gummel scheme. Applying path-following techniques we demonstrate that the model exhibits S-shaped current-voltage curves with regions of negative differential resistance, only recently observed experimentally.