The sublaminate method consists of using stacked and interconnected plates to evaluate interfacial tractions. A high-order plate theory that includes shear and through-thickness stretching is used for each layer. For composites, the stacking sequence information is included. Because the method is an accurate and convenient way to evaluate debond between layers, it is natural to apply the technique to bonded joints. Previous work had focused on exact solutions of these systems. To create a practical tool for bonded joints, nonlinear material properties had to be included. This was accomplished with an approximate method using the P-element technique. One unusual feature is that the material property distribution is approximated using the same functions. The paper outlines the method, and gives examples that highlight the capability of the code. In particular, the bending behavior of joggled joints can be evaluated. The code can also be used to determine strain-energy-release rate for an existing crack between layers.