Surface roughness involved in elastic adhesion of solids may lead to an unstable loading process, which makes experimental characterization and theoretical modeling full of challenge. The full self-consistent model (FSCM) is the most accurate adhesion contact model in the continuum framework, but the strongly nonlinear coupling relationship between surface gap and interaction is very difficult to solve even with high-precision numerical calculation. In this paper, several techniques including the arc-length control method, the Riemann–Stieltjes integral, the adaptive mesh and the Newton–Raphson iterative method were employed in a new solving program to overcome the difficulty of numerical calculation. This calculation method has high efficiency and precision, and it gives a full force–displacement curve containing all possible stable and unstable equilibrium states. The full force–displacement curves for the adhesive contact of wavy surfaces are zigzag and some additional oscillations or minor roundabouts are observed, which have not been detected in the extended JKR and Maugis–Dugdale models. With the increase of surface roughness, the pull-off force first increases and then decreases, because the adhesion mode transforms from the simple contact mode to the multiple one. It is found that the pressure oscillation induced by waviness leads to interface strengthening for small roughness, and the limitation of the theoretical strength and the development of interface cavitation result in interface weakening for large roughness. Compared to the FSCM, the JKR theory significantly overestimates the pull-off forces for large surface roughness. The extended Maugis–Dugdale model is also invalid for large roughness, and if the step cohesive stress is determined by applying a condition of the identical pull-off force at the rigid limit, it can work well for small roughness.