Ion implantation of different dopants (donors and acceptors) into crystalline silicon with subsequent thermal annealing is used for the formation of ultra-shallow p-n junctions in VLSI technology. The experimentally observed phenomenon of transient enhanced diffusion (TED) during annealing hinders further downscaling of advanced VLSI circuits. However, modern mathematical models of dopant diffusion, which are based on the so-called “five-stream” approach, and software packages such as SUPREM4 encounter difficulties in describing TED. In this work, an extended five-stream model for diffusion in silicon is developed, which takes into account all the possible charge states of point defects (vacancies and silicon self-interstitials) and diffusing pairs “dopant atom–vacancy” and “dopant atom–silicon self-interstitial”. The model includes diffusion and drift of differently charged point defects and pairs in the internal electric field and the kinetics of interaction between unlike species. The expressions for diffusion fluxes and sink/source terms that appear in the non-linear, non-steady-state reaction-diffusion equations are derived for both donor and acceptor dopants accounting for multiple charge states of the diffusing species.