The dynamics of flexible fibers or filaments immersed in a fluid are important to understanding many interesting problems arising in biology, engineering, and physics. For most applications, the flows are at very low Reynolds numbers, and the fibers can have aspect ratios of length to radius from a few tens to several thousands. This class of problems is difficult to solve accurately to a reasonable cost with grid based methods, partly due to the different scales in length and radius of the fibers and the fact that elastic equations must be solved within the fibers. Making explicit use both of the fact that we are considering Stokes flow, as well as of the slenderness of the fibers, we have designed a cost-effective method to simulate multiple interacting elastic fibers in a three dimensional Stokes flow. The key points are that for Stokes flow, boundary integral methods can be employed to reduce the three-dimensional dynamics to the dynamics of the two-dimensional fiber surfaces, and that using slender body asymptotics, this can be further reduced to the dynamics of the one-dimensional fiber center-lines. The resulting integral equations include both the effect of the fibers on the flow field, as well as the interactions of fibers, as mediated by the flow. We have developed a numerical method based on this theory that allows for simulating multiple interacting highly flexible fibers. Considering the efficiency of the method, another important fact is that the framework is suitable for introducing a semi-implicit time-stepping scheme, eliminating the severe constraint on the time-step size arising from the elasticity. Our numerical approach is based on second-order divided differences for spatial derivatives, combined with special product integration methods that reflect the nearly singular nature of the integral operators