Abstract Solving three-dimensional space-dependent neutron kinetics equations in hexagonal-z geometry is presented in this paper. The second order backward differentiation formula (BDF) is used to discretize the time derivative terms of the equations. The semi-discretized equations are the form of the fixed source problems, which can be solved by the function expansion nodal method in the hexagonal-z geometry (Cai et al., 2016). The second order accurate time integration BDF with an adaptive step size algorithm is applied. The backward Euler method and the precursor analytic integration method are also implemented and compared with each other. Some numerical tests are used to compare these methods. The numerical results show that BDF is rather more accurate than the backward Euler method and the precursor analytic integration method when the time step is same. The results also show that BDF with adaptive step size is very efficient.