In 1973, Black and Scholes [3] proposed their famous formula for pricing options under the pure-diffusion model. Later Merton [21] proposed to add lognormally distributed jumps, while Kou [18] suggested a model with double exponentially distributed jumps to improve Black and Scholes’ model. In most cases, these models are treated with numerical methods. One of the numerical methods for finding option prices is related to solving a partial integro-differential equation (PIDE). For discretization of this PIDE, most existing methods employ straightforward second-order schemes for spatial direction and time-stepping schemes for time direction. Feng and Linetsky proposed to use the extrapolation approach in combination with implicitexplicit Euler (IMEX-Euler) scheme [9]. Lately, Tangman et al. [27] proposed to use an exponential time integration (ETI) scheme for handling the time direction when solving a PIDE. In Chapter 1, we mainly discuss the history of option pricing problems and numerical methods for pricing options. Amongthem, we mainly focus on Feng and Linetsky’s IMEX extrapolation scheme and Tangman et al.’s ETI scheme. In [19], Lee et al. proporsed a fast approach for computing the Toeplitz matrix [An]j,k = aj−k multiplied by a vector. In Chapter 2, we employs the Toeplitz matrix exponential (TME) method for pricing options. The main idea is using the shift-and-invert Arnoldi method to omit the direct computation of the matrix exponential. In this thesis, we use that method for the option pricing problem. However, the convergence analysis in [19] is not directly applicable in our case. Therefore, we propose another criterion to judge the capability of the shift-and-invert Arnoldi approximation, and prove that such criterion is satisfied in the option pricing problem. Numerical results are given to demonstrate the efficiency of our proposed scheme, with comparison to other numerical methods for the option pricing problem.