En este documento, investigamos un método numérico práctico para resolver una ecuación de difusión fraccional espacial bidimensional unidimensional con coeficientes variables en un dominio finito, basado en el método clásico de Crank-Nicolson (CN) combinado con la extrapolación de Richardson. Se obtienen estimaciones numéricas exactas de segundo orden en el tiempo y el espacio. Se prueba la estabilidad y convergencia incondicional del método. También se presentan y comparan dos ejemplos numéricos con la solución exacta.