We develop a novel real-space Green’s function method for the study of the Cu-substituted Fe-based superconductors based on the Green’s function equation of motion approach. As an extension of the Hubbard-I approximation for the one-band Hubbard model, a series of decoupling scheme are introduced for the multi-orbital Hubbard model. Since the real-space fluctuation induced by the disorder effect can be treated exactly, our method is very suitable for the study of the strongly correlated multi-orbital systems with the presence of impurities. An inhomogeneous three-orbital Hubbard model is employed to involve the electrons in the three dxy , dxz , and dxz orbitals of the iron and cupper atoms in Cu-substituted Fe-based superconductors. When the orbital-selective Mott transition happens in the iron-based superconductors, a Mott gap can be found in the dxy orbital, while the degenerated d xz / yz orbitals remain metallic. With the increase of the Cu substitution, the dopants will induce in-gap states near the Fermi level within the Mott gap of the dxy orbital, and the Mott gap can be fulfilled by the in-gap states at higher doping cases. In addition, we find that the in-gap states only appear in a few discrete lattice points, where the occupancy of dxy orbital is away from half filling. It is also found that all these electronic states are Anderson localized. Meanwhile, the orbital polarization of the whole system increases with the enhancement of the real space fluctuation of the orbital occupancy induced by the Cu substitution. We also find that the Hund’s rule coupling has strong effect in suppressing the real space fluctuation, which causes the in-gap states to vanish. As a result, the dxy orbital will restore to Mott-Hubbard insulator, suggesting that the critical value J c of Hund’s rule coupling for the Mott-Hubbard transition will be pushed up significantly.