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.
