Need to animate contour plots

so I have a code that solves the heat diffusion equation for 2D, I'm pretty sure the function is correct since I had to modify from an existing code for a 1D heat diffusion thing. So I need to animate the contour plot for every nth frame. It's not working out and it looks like the contour plots don't change at all and I can't get the animation to work. Any advice?

import numpy as np
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg')
plt.figure(facecolor='white')


import time
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import animation, rc
from IPython.display import HTML
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg')

# Our FTCS Euler step function will take in the temperature grid 
# array, the x index, and the step sizes for time and space
def ftcseulerstep(T0,i,j,dt,dx):
  T1 = T0[i,j] + c*(T0[i+1,j] + T0[i-1,j] - 2*phi[i,j]) + c*(T0[i,j+1] + T0[i,j-1] - 2*phi[i,j])  
  #row x collum (x,y)
  return T1



# Constants
L = 25     # Mesh grid size
D = 1.9e-5   # Thermal diffusivity constant
N = 25   # Number of divisions in grid
h = (L/N)       # Grid spacing
dt = 10**(-4)     # Time-step size
c = D*dt/(h**2) # Coefficient for euler step function
epsilon = dt/1000
tend = 0.5
t = np.arange(0,tend,dt)
#print(len(t))

#temps#  #row x collum 
Tair = 20
T1 = 50
T2 = 40

phi = np.zeros([N+1,N+1]) + 20 
phi_prime = np.zeros([N+1,N+1]) +20

phi[:,0] = 0
phi[N,:] = 0
phi[0,:] = T1
phi[:,N] = T2

phi_prime[0,:] = T1
phi_prime[:,N] = T2
phi_prime[:,0] = 0
phi_prime[N,:] = 0


phi_prime[:,1] = 0
phi_prime[N-1,:] = 0
phi_prime[:,N-1] = T2
phi_prime[1,:] = T1
#print(phi,phi_prime)


#Store the temperatures 
#grids for every 5000th time step in your solution in a single 3D array.

# dt = 10e-4
# tend = 25
# t = np.arange(0,tend,dt) ---> 25/10e-4 = 250,000
# store at t = 0, 5000, 10,000,... till 250,000  (50 matrixs)


goo =int(( tend/dt ))
doo = (goo/50000)
z = np.zeros([goo,N+1,N+1])  #matrix of row = N+1, collum = N+1 


nsteps = len(t)
for p in range(nsteps):
  for i in range(1,N):
    for j in range(1,N):
      if i == 0 or j == 0 or i == N or j == N:
        phi_prime[i,j] = phi[i,j]
      else:
        #phi_prime[i,j] = ftcseulerstep(phi,i,j,dt,h)
        phi_prime[i,j] = phi[i,j] + c*(phi[i+1,j] + phi[i-1,j] - 2*phi[i,j]) + c*(phi[i,j+1] + phi[i,j-1] - 2*phi[i,j])
        z [p,:,:] = phi
  phi,phi_prime = phi_prime,phi

z_new = z[0:goo:50000]

plt.contourf(phi)

#=====ANIMATION=====#
fig, ax = plt.subplots(figsize=(6,6))
x = np.linspace(0.0,L,N+1)
y = np.linspace(0.0,L,N+1)
X,Y = np.meshgrid(x,y)
Z = X + Y 

ax.set_xlim(( np.min(X), np.max(X)))
ax.set_ylim((np.min(Y), np.max(Y)))
ax.set_xlabel('x [m]',fontsize=16)
ax.set_ylabel('y [m]',fontsize=16)


def animate(n):
    Z_plot = z_new[n] 
    ax.contourf(X,Y,Z_plot,vmin=-1.0,vmax=1.0,cmap='Blues')

    return 

anim = animation.FuncAnimation(fig, animate,
                               frames=(len(t)), interval=40, 
                               blit=False)

# Finally we convert the animation into something we can see
HTML(anim.to_jshtml())