Here we outline an efficient numerical method of solving the time-dependent Schrödinger equation, $i \hbar \frac{\partial | \psi(t) \rangle}{\partial t} = H | \psi(t) \rangle \,.$ The method is based on the Chebyshev polynomial expansion of the time-evolution operator. Only the basic idea is presented in these notes. Further details on the method can be found, for instance, in the following papers:
The Hamiltonian $$H$$ is assumed to have no explicit dependence on time. For instance, if the system under consideration is a particle of mass $$m$$ moving in the presence of a static potential $$V(x)$$, then $\psi(x,t) = \langle x | \psi(t) \rangle$ is the particle's wave function, and $H = -\frac{\hbar^2}{2 m} \frac{\partial^2}{\partial x^2} + V(x) \,.$ In what follows we assume that one can implement the action of $$H$$ on $$\psi$$ numerically. Thus, in the example above, one can discretize the wave function, $\psi(x,t) \quad \to \quad \psi_n(t) \equiv \psi(x_n,t) \qquad (x_n = n \Delta) \,,$ and then approximate the action of $$H$$ by $(H \psi)_n \simeq -\frac{\hbar^2}{2 m} \frac{\psi_{n+1} - 2 \psi_n + \psi_{n-1}}{\Delta^2} + V_n \psi_n \,,$ where $$V_n = V(x_n)$$. Of course, one has to take proper care of the boundary conditions dictated by the system under consideration.
The general solution to the Schrödinger equation can be formally written as $| \psi(t) \rangle = e^{-i H t / \hbar} | \psi(0) \rangle \,,$ where $e^{-i H t / \hbar} \equiv \sum_{n=0}^{\infty} \frac{1}{n!} \left( \frac{t}{i \hbar} \right)^n H^n$ is the time-evolution operator. In practice, however, it is very inefficient to evaluate the operator directly, using the above Taylor expansion. A better approach is to make use of the following expansion: $e^{-i H t / \hbar} = J_0(\tau) I + 2 \sum_{n=1}^{\infty} (-i)^n J_n(\tau) T_n(\mathcal{H}) \,.$ Here, $$I$$ is the identity operator, $$J_n( \cdot )$$ is the Bessel function (of the first kind) of order $$n$$, $$T_n( \cdot )$$ is the Chebyshev polynomial of degree $$n$$, and $$\tau$$ and $$\mathcal{H}$$ are, respectively, the dimensionless time and Hamiltonian, defined as $\tau = \frac{\Vert H \Vert}{\hbar} t \,, \qquad \mathcal{H} = \frac{H}{\Vert H \Vert} \,.$ The norm $$\Vert \cdot \Vert$$ is chosen in such a way that $$\Vert T_n(\mathcal{H}) \Vert \le 1$$ for all $$n$$. As $$|J_n(\tau)| \le |\tau|^n / 2^n n!$$, the series converges rapidly for any fixed $$\tau$$. Thus, truncating the sum at some finite $$n = N$$ (which generally depends on $$\tau$$ and a desired degree of accuracy), one obtains an accurate numerical approximation for the time-evolution operator. The numerical evaluation of the operators $$T_n(\mathcal{H})$$ is straightforward in view of the following recurrence relation: $T_0(\mathcal{H}) = I \,,$ $T_1(\mathcal{H}) = \mathcal{H} \,,$ $T_n(\mathcal{H}) = 2 \mathcal{H} \, T_{n-1}(\mathcal{H}) - T_{n-2}(\mathcal{H}) \,.$
In summary, the numerical approximation of $$| \psi(t) \rangle$$ is given by $| \psi(t) \rangle \simeq J_0(\tau) | \psi_0 \rangle + 2 \sum_{n=1}^N (-i)^n J_n(\tau) | \psi_n \rangle \,,$ where $| \psi_0 \rangle = | \psi(0) \rangle \,,$ $| \psi_1 \rangle = \mathcal{H} | \psi_0 \rangle \,,$ $| \psi_n \rangle = 2 \mathcal{H} | \psi_{n-1} \rangle - | \psi_{n-2} \rangle \,.$