A stress function-based analysis is proposed to provide a simple and efficient approximation method of three-dimensional (3D) state of stress that exists near the free edge of bonded composite patches. In order to apply plane strain assumption in a composite patch, a linear superposition of sliced section from a bonded patch is used. In addition, to describe the load transfer mechanism from the substrate to the composite patch, a simple shear lag model is introduced. The 3D stress behavior at the free edge of the composite patch is modeled by Lekhnitskii stress functions, and the governing equations of the given composite patch are obtained by applying the principle of complementary virtual work. After a suitable expansion of the functions, the governing equations are transformed into two coupled ordinary differential equations, and they are solved by a general eigenvalue solution procedure. As the number of base functions increases, the interlaminar stresses converge. The interlaminar stresses reach maximum at the free edge and decrease sharply at the inner part of the patch. The interlaminar stresses are concentrated at the interface between the layers because of the mismatch of material properties and the geometric singularity. Since the proposed method accurately predicts the 3D stresses in a composite patch bonded on the metal substrate, it can be used as a simple and efficient analytical tool for designing such structural components.