The traditional fixed-point iteration method is typically used for neutronics/thermal-hydralics coupling problems in most nuclear safety analysis codes. But the fixed-point iteration method has a tendency to fail to be used in computing the coupling problems due to slow convergence rate in some cases and even no convergence, and thus resulting in a limited efficiency, especially for the tight-coupling and fast-transient problems. In addition, for the reactor thermal-hydraulic calculation, the traditional finite difference or volume method (FDM or FVM) is used. However, both FDM and FVM require fine mesh size to achieve the desired precision and thus also result in a limited efficiency for the large scale problems. In this paper, to ensure the accuracy, efficiency and convergence for large-scale and complicated coupling problems, the new methods-NEM_JFNK are successfully developed to simultaneously solve the neutronics/thermal-hydralics coupling problems by combining the advantage of the efficient coarse nodal expansion method (NEM) and Jacobian-Free Newton-Krylov method (JFNK). The NEM has been widely used in the reactor physics analysis due to its high efficiency and accuracy in the reactor physics analysis, and it has proved to be superior to FDM and FVM. To improve the efficiency and accuracy for the large scale problems, the NEM is first extended to thermal hydraulic problems from the reactor physics calculation. Then all the governing equations of the neutronics/thermal-hydralics coupling problems can be discretized by the NEM and all the variables can be solved on the coarse meshes so that the size of coupling problems is greatly reduced. To ensure the high accuracy for the coupling problems on the coarse meshes, the high-order coefficients in NEM are successfully transferred between the coupling terms by our research. After that, to ensure the convergence of complicated coupling problems, JFNK based on the NEMs needs to be developed. However, the researches of JFNK based on the NEM in reactor analysis are less and the existing JFNK methods are mostly based on FVM or FDM or the finite element method. In this paper, the NEM discretization equations are successfully integrated into the framework of JFNK through the special treatment and the NEM_JFNK with linear-based preconditioner named LP_NEM_JFNK is also successfully developed. In addition, to take advantage of the existing code and avoid the construction of residual formulations, the non-residual construction NEM_JFNK named NRC_NEM_JFNK is presented and the black-box coupling method is achieved by NRC_NEM_JFNK so that the existing codes only need the simple modification to achieve the combination of the NEM and JFNK. Numerical results show LP_NEM_JFNK and NRC_NEM_JFNK outperform traditional fixed-point iteration method in the sense of convergence rate and efficiency. Further studies are needed to extend the NEM_JFNK method to the multi-dimensional neutronic/thermal hydraulic coupling problems in the high temperature gas-cooled reactor.