Due to environmental concerns, modern transportation solutions demand drastic reductions of fuel consumption and emissions, which can only be achieved with advanced structures using high-performance lightweight materials. For joining these dissimilar materials, adhesive bonding appears as an optimal solution, since mechanical fastening adds weight to the structure, and welding technology is not easily applicable to reinforced plastics and composites. However, one of the major drawbacks associated with bonded joints is the presence of stress concentrations at the overlap ends, especially in single lap joints. In order to reduce these stress concentrations, several techniques have been developed. One of these is the functionally graded adhesive, in which the adhesive properties gradually vary along the overlap length, leading to a more uniform stresses distribution and improving the joint strength. However, the manufacture of an adhesive layer with properties which gradually vary is complex in practice and so is the creation of numerical models that represent these configurations. In the present work, a numerical model for different-graded distribution of adhesive properties along the overlap was developed, using programmed step functions on finite element analysis background in order to discretize and simplify the continuous properties distribution gradient. Cohesive zone modelling was introduced in the numerical model, enabling it to effectively predict graded joint strength. The model was validated with experimental results of functionally graded joints available in the literature. The numerical model developed presents itself as a powerful tool to predict joint strength for functionally graded joints, without imposing large computational demands.