Time integration of Fourier pseudospectral DNS is usually performed using the classical fourth‐order accurate Runge‐Kutta method or other second‐ or third‐order methods, with a fixed step size. We investigate the use of higher‐order Runge‐Kutta pairs and automatic step size control based on local error estimation. We find that the fifth‐order accurate Runge‐Kutta pair of Bogacki and Shampine gives much greater accuracy at a significantly reduced computational cost. Specifically, we demonstrate speedups of 2× to 10× for the same accuracy. Numerical tests (including the Taylor‐Green vortex, Rayleigh‐Taylor instability, and homogeneous isotropic turbulence) confirm the reliability and efficiency of the method. We also show that adaptive time stepping provides a significant computational advantage for some problems (like the development of a Rayleigh‐Taylor instability) without compromising accuracy.