Adhesively bonded composite joints (ABCJs) have been broadly used to connect multimaterials and show their structural and economic advantages compared to traditional bonding methods. However, robust methods are still desired for efficient and accurate lay-wise stress analysis of ABCJs involving multiple boundaries and layers. The purpose of this work was to extend the stress-function variational method for free-edge stress analysis of composite laminates with a finite length. At each interface of the laminate, two unknown Lehknitskii’s stress potential functions were introduced to interpolate the stresses across the layer. A set of 4th-order governing ODEs of the functions was obtained via evoking the complementary virtual work, solved by eigenvalue-function method under proper traction conditions. Corresponding MATLAB™ program was developed and validated by the FEM (ANSYS®). This method can also examine the stress-suppression effect after composite laminates interleafing. Consequently, the above method was furthered for determining the laywise stress distribution in ABCJs.