Two-dimensional time-fractional diffusion equations with given initial condition and homogeneous Dirichlet boundary conditions in a bounded domain are considered. A semidiscrete approximation scheme based on the pseudospectral method to the time fractional diffusion equation leads to a system of ordinary fractional differential equations. To preserve the high accuracy of the spectral approximation, an approach based on the evaluation of the Mittag-Leffler function on matrix arguments is used for the integration along the time variable. Some examples along with numerical experiments illustrate the effectiveness of the proposed approach.