“…u(x,y,t) = 0, t ∈ I, (x,y) ∈ ∂Ω (1.2) u(x,y,0) = ψ 0 (x,y), (x,y) ∈ Ω (1.3) with orders 1/2 < β, γ < 1, constants K x , K y > 0, solution domain Ω = (a,b)×(c,d), and Riesz fractional derivatives Since closed-form analytical solutions of fractional models are rarely accessible in practice, the numerical solutions become very prevalent to empower their successful applications. In literatures, numerical methods of SFDEs proposed to achieve high accuracy and efficiency include finite difference [3][4][5][6][7][8][9][10], finite element (FE) [11][12][13][14][15][16], finite volume [17][18][19] and spectral (element) [20][21][22] methods. It must be emphasized that no matter which discretization is applied, there usually persists intensive computational task in nonlocality caused by fractional differential operators [23].…”