This paper presents a theoretical method to investigate the multiple scattering of electro-elastic waves and dynamic stress around a subsurface cylindrical inclusion in a functionally graded piezoelectric material layer bonded to homogeneous piezoelectric materials. The analytical solutions of wave fields are expressed by employing wave function expansion method, and the expanded mode coefficients are determined by satisfying the boundary conditions around the inclusion. The image method is used to satisfy the mechanical and electrically short conditions at the free surface of the structure. Through the numerical solutions of dynamic stress concentration factors around the inclusion, it is found that when the cylindrical inclusion possesses higher rigidity and greater piezoelectric constant than the two phases of functionally graded materials, the dynamic stress around the inclusion increases greatly. When the distance between the surface of the structure and the inclusion is smaller, the effect of the properties of the inclusion becomes greater. When the cylindrical inclusion possesses lower rigidity and smaller piezoelectric constant than the two phases of functionally graded materials, the maximum dynamic stress shows little difference; however, the variation of the distribution of the dynamic stress around the inclusion is greater. The effect of the properties of the inclusion on the dynamic stress around the inclusion is greater than that on the electric field. The effects of wave frequency on the dynamic stress and electric field are also examined.