Numerical solution of fractional order diffusion problems with homogeneous Dirichlet boundary conditions is investigated on a square domain. An appropriate extension is applied to have a well-posed problem on ${ℝ}^2$ and the solution on the square is regarded as a localization. For the numerical approximation a finite difference method is applied combined with the matrix transformation method. Here the discrete fractional Laplacian is approximated with a matrix power instead of computing the complicated...