We consider an initial-boundary-value problem for a time-fractional diffusion equation with initial condition $u_0(x)$ and homogeneous Dirichlet boundary conditions in a bounded interval $[0, L]$. We study a semidiscrete approximation scheme based on the pseudo-spectral method on Chebyshev–Gauss–Lobatto nodes. In order to preserve the high accuracy of the spectral approximation we use an approach based on the evaluation of the Mittag-Leffler function on matrix arguments for the integration along the time variable. Some examples are presented and numerical experiments illustrate the effectiveness of the proposed approach.