In this paper, an efficient numerical method is proposed to handle two-dimensional fractional diffusion equations on a finite domain. The proposed method combines the product of Legendre wavelet bases for two spatial dimensions and a time direction. The operational matrix of the proposed method is obtained. Tikhonov regularization is employed to stabilize the system in cases where the final linear system of equations is large. The convergence analysis of the method is studied and some numerical examples are presented to investigate the efficiency and accuracy of the method.