El propósito de este documento es desarrollar un esquema numérico para la ecuación de subdifusión fraccional de cuarto orden en dos dimensiones con coeficientes variables y retraso. Utilizando la aproximación de la derivada de Caputo en el tiempo, se logra un método de diferencia finita con precisión de segundo orden en la dirección temporal. La novedad de este documento es introducir un esquema numérico para el problema bajo consideración con coeficientes variables, término fuente no lineal y constante de tiempo de retraso. Los resultados numéricos muestran que las órdenes de convergencia global para las dimensiones espacial y temporal son aproximadamente de cuarto orden en el espacio y de segundo orden en el tiempo.