In this paper we describe a spectrally accurate, unconditionally stable, efficient method using operational matrices to solve numerically two-dimensional advection–diffusion equations on a rectangular domain.The novelty of this paper is to relate for the first time evolution partial differential equations and Sylvester-type equations, avoiding Kronecker tensor products. Furthermore, to reach large times, the calculation of just two matrix exponentials is required, for which we compare different techniques based on Padéʼs approximations, matrix decompositions and Krylov spaces, as well as a new technique which avoids the computation of matrix exponentials. We also illustrate how to take advantage of multiple precision arithmetic.Finally, possible generalizations to non-linear problems and higher-dimensional problems, as well as to unbounded domains, are considered.