We present TimeEvolver, a program for computing time evolution in a generic quantum system. It relies on well-known Krylov subspace techniques to tackle the problem of multiplying the exponential of a large sparse matrix iH, where H is the Hamiltonian, with an initial vector v. The fact that H is Hermitian makes it possible to provide an easily computable bound on the accuracy of the Krylov approximation. Apart from effects of numerical roundoff, the resulting a posteriori error bound is rigorous, which represents a crucial novelty as compared to existing software packages such as Expokit[1]. On a standard notebook, TimeEvolver allows to compute time evolution with adjustable precision in Hilbert spaces of dimension greater than 106. Additionally, we provide routines for deriving the matrix H from a more abstract representation of the Hamiltonian operator. Program summaryProgram Title: TimeEvolverCPC Library link to program files:https://doi.org/10.17632/vvwvng9w36.1Code Ocean capsule:https://codeocean.com/capsule/8431379Developer's repository link:https://github.com/marco-michel/TimeEvolverLicensing provisions: MITProgramming language: C++Supplementary material: An example which demonstrates the computation of time evolution in a concrete physical system.Nature of problem: Computing time evolution in a generic physical quantum system can be reduced to the numerical task of calculating exp(−iHt)v. Here H is the Hamiltonian matrix, which is large and sparse, i corresponds to the imaginary unit, t denotes time and the vector v represents the initial state. A program is needed to perform this computation efficiently. Since the use of approximation methods is unavoidable, it is important to quantify as rigorously as possible the resulting error. Moreover, in order to facilitate the application to various problems in physics, additional functionalities are needed, in particular for forming the Hamiltonian matrix from a more abstract representation of the Hamiltonian operator.Solution method: The program employs known Krylov subspace methods for calculation the exponential of the large sparse matrix (−iHt) times the vector v. The Arnoldi algorithm is used to form the Krylov subspace and exponentiation of the resulting small matrix is achieved by diagonalization. The fact that (−iHt) is anti-Hermitian makes it possible to calculate the error of the Krylov approximation in terms of an easily-computable integral formula. This allows to choose a maximal size of the time step, after which the method is restarted and a new Krylov subspace is formed, while respecting an adjustable error bound. It is rigorous up to inaccuracies of a one-dimensional numerical integral and effects of finite machine precision, for which we also give an estimate. All linear algebra operations are performed with the Intel® Math Kernel Library and Boost is used for numerical integration. The methods for deriving the Hamiltonian matrix rely on a hashtable representation of Hilbert space.