In this study, the authors propose a hybrid computational model for the reliable simulation of electromagnetic‐wave phenomena emerging in graphene structures, which incorporates two variations of an unconditionally‐stable finite‐difference time‐domain algorithm. The new approach features: (a) a dispersive model that relies on the auxiliary differential equation technique and takes into consideration the graphene's surface conductivity, and (b) a modified discretisation algorithm implementing error‐optimised spatial approximations, which is appropriate for graphene‐free mesh nodes. The resulting update equations are devoid of time‐step restrictions, thus formulating a suitable and useful computational framework for the investigation of contemporary graphene‐based problems, which commonly call for densely‐sampled grids and/or prolonged simulations. After theoretically assessing the fundamental features of the proposed algorithm, the propagation properties of surface waves on graphene are extracted numerically and compared with analytical estimations, so that computational validation is provided. In addition, a realistic surface‐wave coupler configuration is simulated and its main characteristics are extracted accurately, exemplifying the method's robustness on more complex setups.