Analyzing the operational states of multiple energy networks (MEN) in multi-energy systems is crucial for ensuring system stability. The dynamic operational characteristics of different energy flows pose challenges for computational analysis. Traditional steady-state methods are inadequate for addressing the dynamics of MEN, especially when dealing with temporal discrepancies between hydraulic and thermal flows in thermal networks (TN) and the heterogeneity between TN and electrical networks. Therefore, this paper proposes a novel holomorphic embedding method (HEM) based on multi-stage decomposition method. The developed HEM constructs a time coefficient matrix and utilize inner-outer loop recursion to handle the time lag between thermal flow and hydraulic flow in the TN. Additionally, we reconstruct a holomorphic matrix, integrating hydraulic flow to bridge thermal and electric power flows, thereby improving the operational heterogeneity among different networks. Real-case simulations show that when the Taylor expansion order in HEM is equal to 4, the proposed method achieves a mere 1% discrepancy from actual operational data, enhancing computational efficiency by 60% compared to the Newton-Raphson method. Moreover, in this real-case, the TN exhibits a maximum delay response time of 180s compared to electrical networks. Exploiting this delay time effectively increases renewable energy generation within multi-energy systems by 961.58kW per day.