An accurate method for time propagation of the coupled Maxwell and time-dependent Kohn-Sham (TDKS) equation is presented. The new approach uses a simultaneous fourth-order Runge-Kutta-based propagation of the vector potential and the Kohn-Sham orbitals. The approach is compared to the conventional fourth-order Taylor propagation and predictor-corrector methods. The calculations show several computational and numerical advantages, including higher computational performance, greater stability, better accuracy, and faster convergence.