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:

An accurate and efficient scheme for propagating the time dependent Schrödinger equation

Unified framework for numerical methods to solve the time-dependent Maxwell equations (in particular, section 3.3)

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 \,.\]