Quantum open systems play an important role in the development of quantum sciences, therefore the study of corresponding numerical method is of great significance. For calculation of quantum open systems, the quasi-adiabatic propagator path integral, which was invented in 1990s, is one of the few numerically exact methods. However, its computational complexity scales exponentially with system size and correlation length, and due to this disadvantage its application is limited in practical calculation. In recent years, the study and application of tensor network has made rapid progress. Representing the path integral by tensor network makes scaling of the computational complexity polynomially, which greatly improve the computational efficiency. Such a new method is called time-evolving matrix product operators. At the very beginning, the reduced density matrix is represented as a matrix product state. Then the time evolution of the system can be achieved by iteratively applying matrix product operators on the matrix product state. The iterative process is amenable to the standard matrix product states compression algorithm, which keeps the computational cost at polynomial scales. The time-evolving matrix product operators is an efficient, numerically exact and fully non-Markovian method, which has a broad application prospect in the study of quantum open systems. For instance, it is already used in the study of the thermalization, heat statistic, heat transfer and optimal control of the quantum open systems, and conversely it can be also used to investigate the effect of the system on the environment. In addition, the TEMPO method is naturally connected to the process tensor which can be used to calculate the correlation function of the system efficiently. This article gives a review of this method and its applications.