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, which is aimed to find an equivalent matrix which is some linearas a linear combinations of Hamiltonians at different time withand the linear coefficient iscoefficients are cleverly chosen to eliminate the error due to Zassenhaus formula. I found this article very useful.