Graphene nanoribbons (GNRs) [1-3] are carbon-based quasi-one-dimensional strips that have provoked great interest nowadays due to some intriguing properties. Simple geometry changes can mainly control their electronic properties. This feature, allied with synthesis advances on tailoring precision, has presented the GNRs as a solid candidate to become a fundamental material in semiconductors applications in the future. Recently, a new edge design approach has gained impulse: heterojunctions. The possibility of mixing completely distinctive edge shapes opened an entirely new fashion in tailoring GNRs. Till now, just a fraction of these materials have their properties fully understood besides the encouraging prospect [2,3]. In this work, we investigated the influence of edge functionalization on arbitrary GNR heterojunctions. The model selected consists of the 2D extended SSH model with electron-phonon coupling [1]. The stationary states rise from the implementation of a self-consistent algorithm, which combines the hamiltonian diagonalization process with the consecutive solving of Euler-Lagrange equations of the expected Lagrangian value. Each heterojunction geometry was subjected to an extensive edge modification. Results unveil a relationship between the stability, electronic properties, and statistical measures of the site's spatial displacement. The energy bandgap tends to be higher when the sites have more energy. Additionally, the more disperse the bond length distribution is, the higher will be the bandgap. Therefore, the connection among these properties has critical importance for future design heterojunctions GNR.