Summary: 3D analytical solutions of stress profiles are obtained for the fibers adjacent to a broken fiber in unidirectional composites with interfacial slip and matrix yielding. To this end, a hexagonal fiber-array model containing a broken fiber is considered in order to derive differential equations based on a shear lag model. By assuming a bilinear stress profile for the broken fiber, and by introducing an elastoplastic shear modulus of the matrix, it is shown that all relevant material parameters are consolidated into a nondimensional characteristic length. The governing differential equations are then analytically solved under the condition that uniform axial deformation prevails in the second and third nearest-neighbor fibers, respectively. The resulting two analytical solutions are verified by numerically solving the governing equations more generally using a finite difference method. The analytical solutions are, moreover, compared with the detailed 3D finite element computations reported recently, leading to the validity of the present solutions and the effectiveness of the nondimensional characteristic length.