Abstract: This paper studies the application of solving high temperature reactor (HTR) neutron diffusion equation with Jacobian-free Newton-Krylov (JFNK) method. Results show that when solving neutron diffusion equation, the relative residual norm of JFNK method decreases slowly at the beginning. Then the rate of convergence become faster and finally reaches a relatively stable value. This feature is conducive to a high-accuracy solution. In the test of two kinds of additional equations, the neutron diffusion equation with flux normalization condition has a better nonlinear convergence behavior. However, due to the longer computational time in solving linear equations, its total computational time is more than the one with k expression. More efficient preconditioning methods should be studied to improve linear equations.