In solving few-group neutron kinetic equations in multidimensions, one must select time step sizes as a function of time such that the temporal truncation error introduced by the discrete time derivative approximation is limited to ensure the desired fidelity. When using the Euler backward finite difference to approximate the first derivative of the flux—a popular approximation because it ensures numerical stability—the truncation error is know to be O(Δt2) and proportional to the second derivative. By employment of the double-time-step-size technique, modified to reduce the frequency that double-time-step-size solutions are required, an estimate of the second derivative can be obtained, leading to an efficient computational algorithm for determining the near-optimum time-step-size sequence to ensure the desired fidelity.