Hi guys I found the answer. The answer is Magnus expansion method. My method is not a good approximation. The method that I described is valid only in a very very small interval of time. Hence the numerical algorithm will be slow and inaccurate. There is a truncation error due to the  Zassenhaus formula 
$$e^{-i/\hbar(t_{1}+t_{2})(H(t_{1})+H(t_{2})}= e^{-i/\hbar\,(t_{1}+t_{2})H(t_{1})}  e^{-i/\hbar\,(t_{1}+t_{2})H(t_{2})} e^{1/\hbar^2\frac{(t_{1}+t_{2})^2}{2} [H(t_{1}),H(t_{2})]} \cdots$$. 
The last additional error factor comes since we deal with time dependent hamiltonian and the Hamiltonian of different time  doesn't commute each other. Hence we have to use the [Magnus expansion][1], which is aimed to find an equivalent matrix as a  linear combinations of Hamiltonians at different time and the linear coefficients are cleverly chosen to eliminate the error due to Zassenhaus formula. I found [this article][2] very useful.


  [1]: http://en.wikipedia.org/wiki/Magnus_expansion
  [2]: http://www.sciencedirect.com/science/article/pii/S0370157308004092