We present a novel theory and computational algorithm for modeling electrical stimulation of nerve fibers in three dimensions. Our approach uses singular perturbation to separate the full 3D boundary value problem into a set of 2D "transverse" problems coupled with a 1D "longitudinal" problem. The resulting asymptotic model contains not one but two activating functions (AF): the longitudinal AF that drives the slow development of the mean transmembrane potential and the transverse AF that drives the rapid polarization of the fiber in the transverse direction. The asymptotic model is implemented for a prototype 3D cylindrical fiber with a passive membrane in an isotropic extracellular region. The validity of this approach is tested by comparing the numerical solution of the asymptotic model to the analytical solutions. The results show that the asymptotic model predicts steady-state transmembrane potential directly under the electrodes with the root mean square error of 0.539 mV, i.e., 1.04% of the maximum transmembrane potential. Thus, this work has created a computationally efficient algorithm that facilitates studies of the complete spatiotemporal dynamics of nerve fibers in three dimensions.