This paper focuses on the application of the Backward Differential Formula (BDF) and Anderson's method for neutron diffusion transient simulations. The B-spline finite element method is used for spatial discretization of continuous equations into the system of ODEs. BDF with adaptive step size and approximation order is used for integration of the ODEs over time. In-step multigroup equation is solved using the block Gauss–Seidel method, where the inner equation is solved using the ILU(0)-preconditioned Conjugate Gradient method. Anderson's method is employed first to solve the in-step non-linear equation that couples the neutron and temperature models, and second, to accelerate the convergence of the block Gauss–Seidel algorithm. The performance of the time integration scheme and the in-step equation solver are verified on the point kinetics equation and two-dimensional benchmarks.