For a new nonlinear iterative method named as Picard-Newton (P-N) iterative method for the solution of the time-dependent reaction-diffusion systems, which arise in non-equilibrium radiation diffusion applications, two time step control methods are investigated and a study of temporal accuracy of a first order time integration is presented. The non-equilibrium radiation diffusion problems with flux limiter are considered, which appends pesky complexity and nonlinearity to the diffusion coefficient. Numerical results are presented to demonstrate that compared with Picard method, for a desired accuracy, significant increase in solution efficiency can be obtained by Picard-Newton method with the suitable time step size selection.