Can I write a Schrodinger equation for time-dependent Hamiltonian like this:
$$i\hbar\frac{d}{dt}\psi(t) = H(t)\psi(t)$$
and then perform Euler integration like this:
$$\psi(t+\Delta t) = (1-\frac{i}{\hbar} H(t)\Delta t)\psi(t)$$
We can do this when $H$ is time-independent for very small $\Delta t$. But when $H$ is time-dependent is it valid to do this?
The Euler integration $$ \psi(t+\Delta t) = (1-\mathrm{i}H\Delta t)\psi(t)\implies \psi(t) = (1-\mathrm{i}H\Delta t)^{t/\Delta t}\psi(0)$$ is not numerically stable even for time-independent or even constant $H$ (that is, an eigenstate) because the normalization of the wavefunction is not preserved - the operator $1-\mathrm{i}H\Delta t$ is not unitary: $$ \lvert \psi(t+\Delta t)\rvert^2 = \lvert (1-\mathrm{i}E\Delta t)\psi\rvert^2 = \lvert \psi\rvert^2 (1+\mathrm{i}E\Delta t)(1-\mathrm{i}E\Delta t) = (1+(E\Delta t)^2)\lvert\psi\rvert^2 \geq \lvert \psi(t)\rvert^2$$ If you try this computationally, your wavefunction will just blow up after enough steps, and it won't be normalized even after a few. That the error is of order $(\Delta t)^2$ will generally not save you unless your $H$ is nicely bounded from above such that you can make $E\Delta t$ sufficently small as well, but this loss of unitarity is really not what you want in any case.
A proper discretization of the Schrödinger equation must use a unitary time-stepping operator, such as $$ U(\Delta t) = \left(1 + \frac{\mathrm{i}}{2}H\Delta t\right)^{-1}\left(1 - \frac{\mathrm{i}}{2}H\Delta t\right)$$ which will work for suitably well-behaved time-independent Hamiltonians (preferably on a spatial lattice as well).
If $H$ is time-dependent, you must choose your time-steps so small such that $H(t+\Delta t)\approx H(t)$, where the necessary accuracy of this approximation depends heavily on your specific Hamiltonian and the chosen numerical approximation method. Nothing more can be said without knowing both of these.
While looking up a ref. for this answer google search produced this thesis on "Numerical methods for solving the TDSE". It is right on the topic, contains all of the following, and more. I am still posting this because it may prove useful as a fast guide.
Stable numerical methods for the TDSE are usually implicit methods. That is, each time-iteration involves numerically solving for a large, but generally sparse, system of linear equations, which means that you need to select adequate solvers. The method mentioned by ACuriousMind is the simplest in this class and is sometimes referred to as the Crank-Nicholson method. It is part of a hierarchy of methods based on Padé approximations to the exponential and involving higher powers of the time step and the Hamiltonian. The accuracy of a diagonal method of order n, using terms up to $(H\Delta t)^n$ on both sides of the equation for $\psi_{n+1}$, is ${\mathcal O}(\Delta t^{2n+1})$.
A cheap alternative, if used carefully, is provided by conditionally stable explicit iteration methods, possibly with adaptive time-step. The price to pay is that the single-step implicit procedure is traded for a multi-step one, involving the solution at two or more steps in time. The simplest such method is known as the second order differencing scheme, see here (paragraph around Eq.(2.175)) and this paper. The idea is to use a 2nd order discretization for the time derivative, which gives $$ \psi_{n+1} - \psi_{n-1} = 2i (H\Delta t) \psi_n $$ The method has the advantage that it conserves the norm (and energy for time-independent Hamiltonians), but stability requires that $$ \Delta t \le \hbar/E_{\text{max}} $$ where $E_{\text{max}}$ is the maximum eigenvalue of the discretized Hamiltonian. For a time-dependent Hamiltonian this means that you need an estimate of $E_{\text{max}}$ over the entire duration of the simulation or some way to adjust the time-step accordingly at each iteration. In addition, the first step requires knowledge not only of $\psi_0$ but also of $\psi_{-1}$, which has to be precalculated by some other suitable method.
More sophisticated methods you may want to consider: Chebyshev iteration, Lanczos reduction.
The short answer: the unitary time-evolution operator in quantum mechanics is
$$ U(0,t) = \hat{T} \exp\left(\frac{i}{\hbar}\int \limits_0^t \mathrm dt' H(t')\right), $$
where $\hat{T}$ denotes time ordering. This is the unitary operator that yields the correct TDSE. The longer answer...
The exponential of a matrix is defined by
$$\exp(\mathbf{M}) = \sum\limits_{n=0}^{\infty} \frac{\mathbf{M}^n}{n!}.$$
(Aside: The exponential of a matrix always converges for finite dimensional matrices.)
Suppose that $\mathbf{M}$ depends on time; i.e. $\mathbf{M} \equiv \mathbf{M}(t)$. Then after a little bit of algebra, one can show that
$$ \frac{\mathrm d}{\mathrm dt}\exp\left(\mathbf{M}(t)\right) = \int\limits_0^1 \mathrm ds \exp\left(s\mathbf{M}(t)\right)\frac{\mathrm d}{\mathrm dt}\mathbf{M}(t)\exp\left((1-s)\mathbf{M}(t)\right). $$
Using Baker-Campbell-Hausdorff,
$$ \frac{\mathrm d}{\mathrm dt}\exp\left(\mathbf{M}(t)\right) = \sum\limits_{n=0}^{\infty} \frac{1}{n!}\left[\mathbf{M}(t),\left[\dots,\left[\mathbf{M}(t),\frac{\mathrm d}{\mathrm dt}\mathbf{M}(t)\right]\right]\dots\right]\exp\left(\mathbf{M}(t)\right), $$
where the $\text{n}^{\text{th}}$ term has $n$ commutators with $\frac{\mathrm d}{\mathrm dt}\mathbf{M}(t)$.
What you want is the following: a unitary time evolution operator $U(0,t)$ that satisfies
$$ \frac{\mathrm d}{\mathrm dt} U(0,t) = \frac{i}{\hbar} H(t) U(0,t), $$
because $U(0,t)\psi(0) = \psi(t)$ for a state $\psi$ (assuming t > 0). Substitution of $\psi(t)$ yields the TDSE. Naively solving the differential equation for $U$ gives $U(0,t) = \exp\left(\frac{i}{\hbar}\int \limits_0^t \mathrm dt' H(t')\right)$. However, expansion of the exponential using the aforementioned formula involves commutators of the Hamiltonian with itself at different times. Because $H$ is an operator, it is not guaranteed to commute with itself at different times. While you can write the time evolution operator in this way, the evolution operator itself will not obey the TDSE, i.e. you will have to use the full matrix expansion to obtain the correct time-dependent evolution equation for the state $\psi(t)$. The correct way to solve this differential equation--such that the unitary operator acting on $\psi$ obeys the TDSE--is recursively. The formal solution is
$$U(0,t) = 1 + \frac{i}{\hbar}\int\limits_0^{t}\mathrm dt' H(t')U(0,t').$$
The differential equation has now been converted to an integral equation. Continue to substitute the solution for $U$ into the above equation. Note that $0 < t' < t$. For the second iteration, you will find that $0 < t'' < t' < t$. This pattern continues ad infinitum. This provides a time-ordering to the expansion, which ensures that the Hamiltonians act on the state $\psi$ in the correct order in time. This yields the time-ordered exponential:
$$ \hat{T} \exp\left(\frac{i}{\hbar}\int \limits_0^t \mathrm dt' H(t')\right) = \sum\limits_{n=0}^{\infty} \left ( \frac{i}{\hbar} \right )^n \int\limits_{t''}^t \mathrm dt' \dots \int\limits_0^{t^{(n-1)}}\mathrm dt^{(n)} H(t')H(t'')\dots H(t^{(n-1)})H(t^{(n)}) $$
It is this unitary operator which a) evolves $\psi$ from an earlier time to a later time and b) obeys the TDSE on its own. Altogether:
$$ \frac{\mathrm d}{\mathrm dt}\psi(t) = \frac{\mathrm d}{\mathrm dt}U(0,t)\psi(0) = \frac{i}{\hbar}H(t)U(0,t)\psi(0) = \frac{i}{\hbar}H(t)\psi(t). $$
NB: the operator without time-ordering is a unitary operator for time-translation, but the time-translation in infinitesimal time-steps zig-zags through time instead of taking a time-ordered path. Because that operator has an expansion in an infinite number of nested commutators, its action on $\psi$ does not yield a TDSE at first glance, but I think it should be equivalent to the action of the time-ordered exponential on the state after an infinite number of algebraic manipulations (correct me if I am wrong on this point). The extra factor of $\frac{1}{n!}$ in the non-ordered exponential (c.f. the formula involving the derivative of $\mathbf{M}(t)$ above) accounts for the over-counting of equivalent paths through time.
The linearized formula for the "evolution by $\Delta t$" is exactly as accurate for a time-dependent $H(t)$ as it is for a time-independent one: the error is of order $O(\Delta t)^2$ in both cases. You must just remember to substitute the right relevant $H(t)$ every time you use the formula at a new time $t$.
It's not hard do see why the accurate doesn't worsen: the extreme change related the time-dependence would be if you substituted $H(t+\Delta t)$ instead of $H(t)$ on the right hand side. But these two again differ by $\Delta t$ times the derivative of $H$ with respect to time (linearization of derivative) and this $H$ and its $O(\Delta t)$ error margin is multiplied by another $\Delta t$ in that equation, so the error in the equation is only $O(\Delta t)^2$ again.
