This paper presents a novel formulation and numerical solutions for adhesively bonded composite joints with non-linear (softening) adhesive behaviour. The presented approach has the capability of choosing arbitrary loadings and boundary conditions. In this model adherends are orthotropic laminates that obey classical lamination theory. The stacking sequences can be either symmetric or asymmetric. Adhesive layer(s) is (are) homogenous and isotropic material. They are modeled as continuously distributed non-linear (softening) tension/compression and shear springs. In this method by employing constitutive, kinematics and equilibrium equations, sets of differential equations for each inside and outside of overlap zones are derived. In the inside of overlap zone, the set of differential equations is non-linear, that is solved numerically. By solving these equations, shear and peel stresses in adhesive layer(s) as well as deflections, stress resultants and moment resultants in the adherends are determined. Most of adhesives have non-linear behavior, therefore unlike previous methods, in which the adhesive layers are modeled as linear materials, in the presented approach the non-linear behavior is assumed for the adhesive layer and can be used to analyze the most of adhesive joints. The numerical results reveal that in the inside of overlap zone, magnitudes of shear forces are considerably large due to high rate of variation in the bending moments. The developed results are successfully compared with those obtained by finite element analysis using ANSYS. The comparisons demonstrate the accuracy and effectiveness of the aforementioned methods.