The coupling-matrix-free iterative s-version finite element method is extended to multiple local meshes in order to model multiple local features such as holes, inclusions, and cracks. The formulation with multiple local meshes is presented, along with the stress transfer method between the local meshes. The present method does not require the generation of coupling stiffness matrices with a very sophisticated numerical integration method between the global and local meshes or between the local meshes. Instead, stress transfers between pairs of overlapping meshes are performed. Several numerical examples, such as a patch test, a two-hole problem, and a structure with many holes, are presented. The examples demonstrate that the present method is capable of representing a uniform stress distribution as well as capturing the stress concentration of multiple holes that are located in the vicinity of each other. Moreover, numerical integration methods of the global mesh and local meshes are found to have significant influences on the convergence of the iteration. In these numerical examples, the straightforward Gaussian quadrature could not achieve convergence, whereas the element subdivision technique could. Thus, sufficient element subdivisions in the global mesh and the local meshes are necessary in order to produce a converged solution.