Summary: In this paper, a time-space fractional diffusion equation in two dimensions (TSFDE-2D) with homogeneous Dirichlet boundary conditions is considered. The TSFDE-2D is obtained from the standard diffusion equation by replacing the first-order time derivative with the Caputo fractional derivative

${}_{t}{D}_{*}^{\gamma}$,

$\gamma \in (0,1)$, and the second-order space derivatives with the fractional Laplacian

$-{(-{\Delta})}^{\alpha /2}$,

$\alpha \in (1,2]$. Using the matrix transfer technique proposed by

*M. Ilić* et al. [Fract. Calc. Appl. Anal. 9, No. 4, 333–349 (2006;

Zbl 1132.35507)], the TSFDE-2D is transformed into a time fractional differential system as

${}_{t}{D}_{*}^{\gamma}\mathbf{u}=-{K}_{\alpha}{\mathbf{A}}^{\alpha /2}\mathbf{u}$, where

$\mathbf{A}$ is the approximate matrix representation of

$(-{\Delta})$. Traditional approximation of

${\mathbf{A}}^{\alpha /2}$ requires diagonalization of

$\mathbf{A}$, which is very time-consuming for large sparse matrices. The novelty of our proposed numerical schemes is that, using either the finite difference method or the Laplace transform to handle the Caputo time fractional derivative, the solution of the TSFDE-2D is written in terms of a matrix function vector product

$f\left(\mathbf{A}\right)\mathbf{b}$ at each time step, where

$\mathbf{b}$ is a suitably defined vector. Depending on the method used to generate the matrix

$\mathbf{A}$, the product

$f\left(\mathbf{A}\right)\mathbf{b}$ can be approximated using either the preconditioned Lanczos method when

$\mathbf{A}$ is symmetric or the

$\mathbf{M}$-Lanzcos method when

$\mathbf{A}$ is nonsymmetric, which are powerful techniques for solving large linear systems. We give error bounds for the new methods and illustrate their roles in solving the TSFDE-2D. We also derive the analytical solution of the TSFDE-2D in terms of the Mittag-Leffler function. Finally, numerical results are presented to verify the proposed numerical solution strategies.