Integration factor (IF) methods are a class of efficient time discretization methods for solving stiff problems via evaluation of an exponential function of the corresponding matrix for the stiff operator. The computational challenge in applying the methods for partial differential equations (PDEs) on high spatial dimensions (multidimensional PDEs) is how to deal with the matrix exponential for very large matrices. Compact integration factor methods developed in [Nie et al., Journal of Computational Physics, 227 (2008) 5238-5255] provide an approach to reduce the cost prohibitive large matrix exponentials for linear diffusion operators with constant diffusion coefficients in high spatial dimensions to a series of much smaller one dimensional computations. This approach is further developed in [Wang et al., Journal of Computational Physics, 258 (2014) 585-600] to deal with more complicated high dimensional reaction-diffusion equations with cross-derivatives in diffusion operators. Another approach is to use Krylov subspace approximations to efficiently calculate large matrix exponentials. In [Chen and Zhang, Journal of Computational Physics, 230 (2011) 4336-4352], Krylov subspace approximation is directly applied to the implicit integration factor (IIF) methods for solving high dimensional reaction-diffusion problems. Recently the method is combined with weighted essentially non-oscillatory (WENO) schemes in [Jiang and Zhang, Journal of Computational Physics, 253 (2013) 368-388] to efficiently solve semilinear and fully nonlinear convection-reaction-diffusion equations. A natural question that arises is how these two approaches may perform differently for various types of problems. In the first part of this dissertation, we study the computational power of Krylov IF-WENO methods for solving high spatial dimension convection-diffusion PDE problems (up to four spatial dimensions). Systematical numerical comparison and complexity analysis are carried out for the computational efficiency of the two different approaches. We show that although the Krylov IF-WENO methods have linear computational complexity, both the compact IF method and the Krylov IF method have their own advantages for different type of problems. This study provides certain guidance for using IF-WENO methods to solve general high spatial dimension convection-diffusion problems.In the second part of this dissertation, we combine the Krylov integration factor methods with sparse grid combination techniques and solve high spatial dimension convection-diffusion equations such as Fokker-Planck equations on sparse grids. Numerical examples are presented to show that significant computational times are saved by applying the Krylov integration factor methods on sparse grids.