Reduced differential transform method

hi , i want to solve in python this equation :$$partial^2_{xx}u+partial^2_{yy}u=partial^2_{xx}u_t-u(1-u) where $(x,y,t)inOmega timesmathbb{R_{+}}$ where $Omega=(-1,1)times (-1,1)$such that $u(x,y,0)=exp(-{x+y))$ such that $u=0 inpartialOmega$
so to do that i use reduced diffrential transform method :
(k+1)U_{k+1)(x,y)=partial^2_{xx}U_{k}(x,y)+partial^2_{yy}U_{k}+U_k(x,y)-sum_{r=0}^{k}U_{r}(x,y)U_{k-r}(x,y)$

and by solving that for $k =0;dot N$ for N large for example here i will take only 10
then $u^{tilde}=sum_{k=0}^{k=N}U_{k}(x,y)t^{k}$ and this is the approximate solution ,i use this code :

import numpy as np
x_min, x_max = -1, 1
y_min, y_max = -1, 1
t_max = 1
N = 10
def initial_condition(x, y):
    return np.exp(-(x + y))
def solve_equation(N, t_max):
U = [initial_condition]
    for k in range(1, N+1):
        def U_k(x, y):
            return np.gradient(np.gradient(U[k-1](x, y), axis=0), axis=0) + 
                   np.gradient(np.gradient(U[k-1](x, y), axis=1), axis=1) + 
                   U[k-1](x, y) - 
                   np.sum(U[r](x, y) * U[k-r](x, y) for r in range(k))
        U.append(U_k)
def u_tilde(x, y, t):
    return np.sum(U[k](x, y) * t**k for k in range(N+1))
return u_tilde
 Compute the approximate solution for given x, y, and t
u_tilde = solve_equation(N=10, t_max=1)
x = np.linspace(x_min, x_max, 100)
y = np.linspace(y_min, y_max, 100)
X, Y = np.meshgrid(x, y)
T = 0.5  # Example time
Z = u_tilde(X, Y, T)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('u(x, y, t)')
plt.show()

but it doesnt work?