A mathematical procedure is developed to utilize the complementary energy method, by minimization, in order to obtain an approximate analytical solution to the 3D stress distributions in bonded interfaces of dissimilar materials. The stress solutions obtained predict the stress jump, at the interfaces, which cannot be captured by current FEA methods. As a novel method, the penalty function is used to enforce the displacement boundary conditions at the interfaces. Furthermore, the mathematical procedure developed enables the integration of different interfacial topographies into the solution procedure. In order to incorporate the effects of surface topography, the interface is expressed as a general surface in Cartesian coordinates, i.e. F(x, y, z) = 0. In this paper, the flat interface problem, i.e. y = 0 surface is considered for verification of the method by comparison with the FEA method. A comparison of the results reveals our new mathematical procedure to be a promising and efficient method for optimizing interface topographies.