The field of option pricing has grown very rapidly over the last 30 years. Many academic researchers and financial practitioners in disciplines such as mathematics, finance, economics, business, physics, computer, etc., have been working on this interesting area. In this PhD thesis, the author, from a mathematics discipline, employs modern Krylov subspace methods with preconditioning techniques to solve a family of Toeplitz-like systems from option pricing. The thesis is divided into three parts. First, we introduce the definition of a family of generating functions (FGF) [119] of Toeplitz matrices, which is a generalization of the generating function of the Toeplitz matrix. The FGF has an important application in pricing derivatives. The pricing of a European call option under Merton’s jump-diffusion model, under certain assumptions, can be described by a partial integro-differential equation (PIDE) without a convection term. We then use the FGF to analyze the convergence rate of the preconditioned conjugate gradient (PCG) method with Strang’s circulant preconditioner which is useful for solving Toeplitz systems arising from a discretization of the PIDE. We show that if the FGF meets certain conditions, then the spectrum of the preconditioned matrix is clustered around 1, and the smallest eigenvalue of the preconditioned matrix is uniformly bounded away from 0. It follows that the convergence rate of the PCG method is superlinear. The FGF provides a rigorous theoretical basis for the convergence analysis of the conjugate gradient (CG) method. All aspects related to the generating function could be extended to the FGF. Second, we solve a nonsymmetric Toeplitz system resulting from the finite difference discretization of a PIDE, which is led by a European call option under Merton’s jump-diffusion model. In Merton’s model, the jump magnitude distribution is normal with mean µJ and standard deviation σJ . When µJ = 0, discretizing the PIDE without the convection term yields a symmetric Toeplitz system [102, 119]. While for µJ 6= 0, the resulting system is a nonsymmetric Toeplitz system. In [102, 119], the authors only consider the case µJ = 0. In this part, we discuss a more general case with µJ 6= 0. The conjugate gradient normal equation residual method [53, 101] is employed. To speed up the numerical process, preconditioning techniques are used. We apply the CG method to a normalized preconditioned system [29]. Strang’s circulant preconditioner [111] and a tri-diagonal preconditioner are studied. By using the definition of the FGF introduced in [119], we show that the spectrum of the normalized preconditioned matrix with Strang’s preconditioner is clustered around 1 and the smallest eigenvalue is uniformly bounded away from 0. It follows that the convergence rate of the CG method when applied to such a normalized preconditioned system is superlinear [25]. Similar results are also obtained for the tri-diagonal preconditioner. We see from the numerical results in Section 6.2 that both preconditioners work very well. Third, we exploit quadratic finite element (FE) with preconditioning for option pricing problems in the Stochastic Volatility with correlated and Contemporaneous Jumps in return and variance (SVCJ) model. The value function for a European vanilla option or barrier option satisfies a PIDE. After spatial discretizations of a variational formulation of the PIDE by quadratic FE, we obtain an ordinary differential equation (ODE) system (which is usually referred to as a semi-discretization of the variational problem). An implicit-explicit Euler based extrapolation scheme proposed by Feng and Linetsky [51] is used to integrate this ODE system in time. For the solutions of the resulting linear systems we adopt the bi-conjugate gradient stabilized (BICGSTAB) method. Since the coefficient matrix of the resulting linear system is block penta-diagonal with penta-diagonal blocks and has certain block structure, a Schur complement preconditioner is designed to accelerate the convergence of the BICGSTAB method. Numerical experiments demonstrate that the combination of quadratic FE, extrapolation scheme, and preconditioning techniques is very efficient to the option pricing problems in the SVCJ model. We expect this combination to be very attractive to financial engineering modelers. In summary, the major contributions of the thesis are: the introduction of the definition of FGF for Hermitian Toeplitz systems, and the FGF provides a rigorous theoretical basis for the convergence analysis of the CG method; the generalization of the definition of FGF to non-Hermitian Toeplitz systems, and the FGF provides a rigorous theoretical basis for the convergence analysis of the CG method when applied to a normalized preconditioned system, which is led by a European call option under Merton’s jump-diffusion model with µJ 6= 0; the use of the quadratic FE method and the Schur complement preconditioner for European option pricing problems under the SVCJ model.