Modeling elastic fields in materials with inclusions is pivotal to understanding the mechanical behavior of composites, yet conventional finite element formulations often suffer from high computational cost when dealing with complex inclusion geometries. In this talk, we present a hybrid Trefftz finite element method based on a single-functional formulation to deal with this problem. The method assumes two independent displacement fields, one within the element domain and the other along its boundary, and couples them through a modified functional with the matrix–inclusion continuity conditions inherently satisfied by the intra-element field. The single hybrid functional couples the intra-element and boundary fields through boundary-only integration, forming a concise and efficient numerical scheme that requires no integration along material interfaces. Numerical studies are carried out on randomly distributed inclusions or voids, to confirm that the proposed method captures inclusion-induced local stress concentrations with high fidelity even on coarse meshes, and significantly reduces degrees of freedom and computational time compared with conventional finite element method. The proposed method can be further applied to problems of multilayer-coated fibers, demonstrating its robustness and efficiency in dealing with complex inclusion configurations. This study provides a new pathway to solve the elastic field of composite materials with high accuracy and high efficiency.