Because charge carriers of many organic semiconductors (OSCs) exhibit fractional drift diffusion (Fr-DD) transport properties, the need to develop a Fr-DD model solver becomes more apparent. However, the current research on solving the governing equations of the Fr-DD model is practically nonexistent. In this paper, an iterative solver with high precision is developed to solve both the transient and steady-state Fr-DD model for organic semiconductor devices. The Fr-DD model is composed of two fractional-order carriers (i.e., electrons and holes) continuity equations coupled with Poisson’s equation. By treating the current density as constants within each pair of consecutive grid nodes, a linear Caputo’s fractional-order ordinary differential equation (FrODE) can be produced, and its analytic solution gives an approximation to the carrier concentration. The convergence of the solver is guaranteed by implementing a successive over-relaxation (SOR) mechanism on each loop of Gummel’s iteration. Based on our derivations, it can be shown that the Scharfetter–Gummel discretization method is essentially a special case of our scheme. In addition, the consistency and convergence of the two core algorithms are proved, with three numerical examples designed to demonstrate the accuracy and computational performance of this solver. Finally, we validate the Fr-DD model for a steady-state organic field effect transistor (OFET) by fitting the simulated transconductance and output curves to the experimental data.