I'm using 4th order RK to solve the schroedinger equation in atomic units. Say I want to simulate 400fs in intervals of h=10fs, then in atomic units this is h=413a.u and 400fs=16500a.u. 4RK involves repeated multiplication by the timestep, so with h=413 everything blows up. Could someone explain what I'm missing?
Edit: to be clear, I'm solving
\begin{align} \frac{d}{dt} \psi = -i H \psi \end{align} Where $\psi$ is a vector and $H$ a hermitian matrix.
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).
