A new three-dimensional (3D) interface element has been developed to study the nature and properties of the resistance to fiber-matrix interfacial debonding in composites composed of ductile fibers in a brittle or elastic matrix. In order to model various types of interface characteristics, an experimental/assumed nonlinear bond stress-slip relationship is used to calculate the stiffness of the discrete pull out springs connecting the steel fiber and the concrete element nodes. The cementitious matrix and the steel fiber are discretized by using 3D eight-noded isoparametric hexahedron elements and two-noded bar elements, respectively, with three translational degrees of freedom at each node. Other 3D elements such as 20-noded quadratic hexahedron, four-noded tetrahedron, six-noded linear, and 15-noded quadratic prism elements, and three-noded quadratic bar elements are incorporated in the program. The von Mises and Mohr-Coulomb yield criteria are used in the material nonlinearity. An incremental total Lagrangian formulation for concrete and an updated Lagrangian formulation for steel fiber and the interface element, are employed in the geometric nonlinearity with large deformation. The effects of fiber size and shape such as straight, hooked-end, and crimped fibers, and the influence of number of fibers being pulled out simultaneously from the same area are examined. Apparently, this study constitutes the first in which the fiber pull-out behavior in fiber-reinforced composites has been studied by using a novel approach in the nonlinear 3D interface element.