Q:

Normalizing a density matrix at each iteration

I need to numerically evolve a density matrix using this formula(Actually I have more terms but right nows I am starting with this and facing problems):
$$\dot\rho(t) = -i[H(t), \rho(t)]$$
$H(t)$ is the Hamiltonian which is of a form $H_o + \epsilon(t)H_1$ where both $H_o$ and $H_1$ are constant $20\times20$ matrices and $\epsilon(t)$ is a piece-wise constant function.  
So, my numerical scheme is basically
$$\rho(t + dt) = \rho(t) -i[H(t), \rho(t)]dt$$
and substituting $\rho(t + dt)$ to the right hand side to get $\rho(t + 2dt)$ and so on.
The problem I am facing is that the entries of the density matrix quickly blow up. 
I have tried checking whether the density matrix remains a valid density matrix at each iteration: namely if it is Hermitian, if it's trace is one and if it's eigenvalues sum to one. 
Now the first two conditions are always met but it's eignevalues are becoming larger and larger. In some sense it seems reasonable too as I could have changed the time-step $dt$ to say 2 times it's value and the eigenvalues of $\rho(t+dt)$ would change. So is it necessary to employ some sort of normalization procedure to the density matrix after each iteration?
Note: I have tried to use very small values of $dt$ and there is no improvement so the problem is not in the value of the time-step.

A:

As always with numerical problems, it's a good idea to simplify things as much as possible to isolate what's really going wrong. You can simplify the problem by considering the evolution of a pure state $\psi$ instead of an entire density matrix, so solving
$$i\dot\psi = H\psi$$
for the vector $\psi$ instead of the matrix $\rho$. Additionally, you can forget the perturbation $\epsilon$ and just assume that $H = H_0$.
The numerical method you're using is called the explicit Euler scheme. The explicit Euler scheme is guaranteed to give a sequence of solutions that, for any instant in time, will converge to the true solution in the limit as $\delta t$ goes to 0. However, it is not guaranteed to preserve every property that might be useful to you, for example the fact that the probability amplitude should be conserved.
To see this in more detail, let's compare the exact solution $\Psi$ with the explicit Euler solution $\psi$. At time $t + \delta t$, you can write down an expression for $\Psi$ using the matrix exponential:
$$\Psi(t + \delta t) = e^{-i\delta t\cdot H}\Psi(t).$$
Since the Hamiltonian $H$ is Hermitian, the operator $e^{-i\delta t\cdot H}$ is unitary. One consequence is that the norm of $\Psi(t + \delta t)$ is the same as the norm of $\Psi(t)$, so the probability amplitude is conserved for the exact system. With the explicit Euler approximation, however, we can write the value $\psi_{n + 1}$ of the wave function at time $(n + 1)\delta t$ in terms of the previous value $\psi_n$ as:
$$\psi_{n + 1} = \psi_n - i\delta t H \psi_n = (I - i\delta t\cdot H)\psi_n.$$
The operator $I - i\delta t\cdot H$ is not unitary, so you can only expect conservation of probability in the limit as $\delta t \to 0$ and $n \to \infty$ such that $n\delta t$ stays fixed. This is really unfortunate because our approximation scheme doesn't preserve one of the most important properties of the physical system.
Luckily, there are many more timestepping schemes besides explicit Euler. Another approach might be to approximate $e^{-i\delta tH}$ by the 1, 1-Pade approximant:
$$e^{-i\delta tH} \approx \left(1 + \frac{i\delta tH}{2}\right)^{-1}\left(1 - \frac{i\delta tH}{2}\right).$$
You can show with a little linear algebra that this approximation is a unitary operator, in which case all of your discrete approximations to $\Psi$ will conserve probability up to floating-point error. This approximation is equivalent to using the Crank-Nicholson discretization of the problem instead of explicit Euler.
Leveque's book is a good introductory reference on finite difference methods. This will give you a lot of the necessary background on constructing stable discretizations of ODE and PDE. It's possible to design numerical methods explicitly for preserving things like the total probability amplitude; these are called symplectic methods. The best reference I know is Hairer, Lubich and Wanner, but unfortunately this book is quite advanced. You could start by looking at the Verlet method, which is the simplest symplectic method for Hamiltonian systems.

