We describe a new and efficient method for numerical study of the dynamics and statistics of single-electron systems presenting arbitrary combinations of small tunnel junctions, capacitances, and voltage sources. The method is based on the numerical solution of a master equation describing the time evolution of the probabilities of the electric charge states of the system, with iterative refining of the operational set of states. The method is able to describe very small deviations from the “classical” behavior of a system, due to finite speed of applied signals, thermal activation, and macroscopic quantum tunneling of charge (cotunneling). As an illustration, we briefly study the leakage rate in single-electron traps and the accuracy of several devices (turnstile, pump, and a hybrid circuit) suitable as standards of dc current.