Abstract: With increasing uncertainties in the source and load, uncertain energy flow calculation is important for optimal energy scheduling of the integrated electricity and natural gas system. A holomorphic embedding affine energy flow method for natural gas network is proposed to address the issues of the complexity of inverting the affine matrix and difficulties in applying it to the looped network of the affine power flow method for natural gas network. On this basis, the energy flow equations of the electricity and gas subsystems are restructured according to the energy flow at the coupling devices, and a holomorphic embedding electricity-gas affine multi-energy flow model is constructed. An affine energy flow alternating recursion calculation method is proposed to reduce the dependency on initial values and improve the efficiency of multi-energy flow calculation. The germ solution of the status variables of the no-load conditions of the power grid is used as the starting point for multi-energy flow calculation. Based on the recursion relationships of the power series coefficients of affine status variables, the power series coefficients of affine status variables of the power grid, gas network, and coupling devices are calculated alternately and recursively, and the energy flow distribution interval of integrated electricity and natural gas system is solved. Simulation results verify the advantages of the proposed method in terms of conservatism, convergence, and computational efficiency.