We present an approach for numerically simulating the dynamics of flexible fibers in a three-dimensional shear flow using a scalable immersed boundary (IB) algorithm based on Guermond and Minev's pseudo-compressible fluid solver. The fibers are treated as one-dimensional neutrally-buoyant Kirchhoff rods that resist stretching, bending, and twisting, within the generalized IB framework. We perform a careful numerical comparison against experiments on single fibers performed by S.G. Mason and co-workers, who categorized the fiber dynamics into several distinct orbit classes. We show that the orbit class may be determined using a single dimensionless parameter for low Reynolds flows. Lastly, we simulate dilute suspensions containing up to hundreds of fibers using a distributed-memory computer cluster. These simulations serve as a stepping stone for studying more complex suspension dynamics involving aggregation of fibers (or flocculation) and particle sedimentation due to added mass. (C) 2015 Elsevier B.V. All rights reserved.