Abstract: The fault propagation in integrated electricity-gas system(IEGS) exhibits complex features of multi-scale and bidirectional dynamics. Traditional time-domain simulation is difficult to achieve both stable and fast cross-scale numerical analysis, resulting in extremely low efficiency of IEGS stability analysis. To conquer this barrier, a variational neural ordinary differential equation (V-NODE)-based long-term dynamic stability inference method for post-faulted IEGS is proposed. First, involving the interactive responses between power grid and natural gas network, a full-system dynamic simulation model for IEGS is established. Then, synthetic minority over-sampling technique is conducted to settle post-fault response dataset with balanced stable and instable trajectories to prevent NODE overfitting in the unbalanced sample space. Furthermore, variational autoencoder is utilized to embed the steady-state variables of the power grid into time-domain. In this sense, inadequate representability and weak generalizability of conventional NODE on multiple steady equilibrium points can be circumvented. Finally, a V-NODE cost-sensitive learning method adapted to multi-stability modes is proposed to prevent NODE from overfitting the out-of-sync trajectories. The improved IEGS case shows that, compared to traditional ODE-based simulator, our method improves efficiency beyond 3 orders of magnitude, with a second elapsed time. It also beats the data-driven rival with substantial superiority in accuracy. Besides, the reachability analysis considering the V-NODE prediction error has verified the effectiveness and trajectory extrapolation ability of the proposed method.