Abstract: High-temperature gas-cooled reactors (HTGRs), as one of the recommended reactor types for Generation Ⅳ advanced nuclear energy systems, are characterized by high thermal efficiency and inherent safety, thus having great application prospects. However, water ingress and air ingress accidents in HTGRs are accidents that require special consideration, in which air or water vapor reacts with the graphite of the reactor core, generating a variety of gases and corroding the graphite. The corrosion of the graphite of fuel pebbles may lead to leakage of radioactive material, and the corrosion of the core components may lead to a reduction in the mechanical strength of the core. To research about the chemical corrosion phenomenon, the research of flow and diffusion process of multiple gases in the core is a prerequisite. And the latter is a more complex nonlinear coupled system than the single-component problem which needs to be solved stably, accurately and efficiently. Newton-Krylov (NK) method has excellent computational performance because it solves simultaneously all physical quantities. This method is more stable, accurate and efficient than Picard method or operator splitting method in solving large nonlinear systems. Its applications to HTGRs include neutronic calculations, thermal-hydraulic calculations and even coupled calculations. And it has been applied to the core, one-loop and two-loop of HTGRs. In this paper, a multi-component gas flow coupling calculation program for the characteristics of HTGRs was developed based on NK method. Among the two variants of the NK method, finite-Difference Jacobian Newton-Krylov method (DJNK) which explicitly constructs the Jacobian was chosen as the solution method. And a matrix-coloring algorithm was used to significantly reduce the number of finite difference times and improve the computational efficiency of DJNK method. Besides, a modified NK method based on the physical properties is proposed to avoid the intermediate iterative values of the non-physical solutions. By studying the selection of preconditioning matrix structure and the number of filling levels of ILU decomposition, the preconditioning matrix structure of “main coupling type” is proposed and ILU(5) is selected as suitable preconditioning for this research, which improves the efficiency of NK method for solving linear equations. The correctness of this program was verified by comparison with the high-temperature gas-cooled reactor analysis program TINTE. The maximum error between this program and TINTE is less than 0.1%. The test results show that the computational performance of NK method is about 6-7 times that of Picard method.