High concentrations of interfacial stress near the adherend ends are primarily responsible for the debonding failure of bonded joints, such as those structured extensively in civil and structural engineering; aeronautical, ground, and marine vehicles; and flexible electronics and microelectronic packaging. Accurate determination of these interfacial stresses is crucial to improved structural design and optimization as well as health monitoring of the structures in which such joints are found. A variety of joint models have been available in the literature for joint strength analysis and structural design. Yet, a few deficiencies still exist in most of these models in accurate prediction of joint stresses, including the violation of the generalized Hooke's law of the adhesive layers and failure to satisfy the physical traction conditions at the free edges of joint adherends. In this thesis, a generalized stress-function variational method is developed for the determination of the interfacial shear and normal stresses in general bonded bimaterial joints subjected to mechanical and thermomechanical loads. Specifically, three types of joints are considered in this study, including single-sided bonded joints, single-sided strap joints, and single-lap joints. During the formulation, two unknown interfacial stress functions are directly introduced to satisfy the traction boundary conditions of the joints; the Euler-Bernoulli elementary beam theory and 2D elasticity are used to determine the stress components of the adherends in terms of the interfacial stress functions. By utilizing the theorem of minimum complementary strain energy, the governing equations of the bimaterial joint are obtained as a system of two coupled 4th-order ordinary differential equations (ODEs) of the introduced stress functions. These ODEs are formatted into a generalized eigenvalue problem, and are further solved numerically by designing robust and efficient computational codes using MATLAB™. The results of the analysis are validated by comparison with elementary mechanics of materials as well as detailed finite element analysis (FEA) using ANSYS™; the current models can accurately satisfy the shear stress-free boundary conditions at the adherend edges. In addition, the proposed method is further applied to the analysis of progressive cracking in hard coatings. In this analysis, a cracked hard coating layer bonded onto a substrate is modeled as a single-sided bonded joint, and the expressions of strain energy derived in this study are incorporated into energy-based cracking criteria of the system. Using the above variational method, the critical loads (i.e., applied axial stress, shear force, bending moment, or temperature change) for initial cracking can be determined. Thus, the present method is also capable of modeling progressive cracking in hard coating systems.