Runge-Kutta timestep in atomic units

This phenomena of an explicit Runge-Kutta method is known as stability. For a given equation, the eigenvalues of your system (here, the Hamiltonian) determine a maximum time step for which the numerical integration scheme will be stable (i.e. not blow up). Here's a good resource which mentions the stability region for RK4. As these resources show, there are special methods such as implicit methods which can be used in order to take larger time steps.

But even explicit Runge-Kutta methods can be enhanced. For example, using adaptive Runge-Kutta methods with step-stabilization (PI-adaptive step controllers) will allow the method to only take smaller time steps when they are needed for stability, and this then allows the steps to be (usually) larger on average. Standard ODE solver software like those in MATLAB, Fortran, or Julia make use of this. So I would recommend using whatever ODE solvers you have available in your language, since not only will they have (potentially PI-controlled) adaptive methods, but also integrators for stiff systems. These will be much more efficient than RK4 which is really a method just for teaching and not for production use (other explicit RK methods have larger stability regions and/or less error for the same amount of work).