2D Heat Conduction with Python

the output of the study///// I just started learning Python, so I am new with python. I have written a simple code for 2D Heat Conduction. I don't know what is the problem with my code. The result is so strange.I think the Temperature distribution is not shown correctly. I have searched about it a lot but unfortunately I could not find any answer for my problem . can anyone help me? Thank you

# Library

import numpy                           
from matplotlib import pyplot     



# Grid Generation

nx = 200
ny = 200                                    
dx = 2 / (nx-1)
dy = 2 / (ny-1)

# Time Step

nt = 50                                                                   
alpha = 1                                 
dt = 0.001                     

# Initial Condition (I.C) and Boundry Condition (B.C)

T = numpy.ones((nx, ny))                         # I.C (U = Velocity)
x = numpy.linspace(0,2,nx)                       # B.C
y = numpy.linspace(0,2,ny)                       # B.C

Tn = numpy.empty_like(T)                         #initialize a temporary array
X, Y = numpy.meshgrid(x,y)

T[0, :] = 20          #  B.C
T[-1,:] = -100        #  B.C
T[:, 0] = 150         #  B.C
T[:,-1] = 100         #  B.C
# Solver
###Run through nt timesteps
    
for n in range(nt + 1): 
    Tn = T.copy()
        
    T[1:-1, 1:-1] = (Tn[1:-1,1:-1] + 
                        ((alpha * dt / dx**2) * 
                        (Tn[1:-1, 2:] - 2 * Tn[1:-1, 1:-1] + Tn[1:-1, 0:-2])) +
                        ((alpha * dt / dy**2) * 
                        (Tn[2:,1: -1] - 2 * Tn[1:-1, 1:-1] + Tn[0:-2, 1:-1])))
        
    T[0, :] = 20          # From B.C
    T[-1,:] = -100        # From B.C
    T[:, 0] = 150         # From B.C
    T[:,-1] = 100         # From B.C

   
fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, T)
pyplot.colorbar()
pyplot.contour(X, Y, T)
pyplot.xlabel('X')
pyplot.ylabel('Y');

    

You are using a Forward Time Centered Space discretisation scheme to solve your heat equation which is stable if and only if alpha*dt/dx**2 + alpha*dt/dy**2 < 0.5. With your values for dt, dx, dy, and alpha you get

alpha*dt/dx**2 + alpha*dt/dy**2 = 19.8 > 0.5

Which means your numerical solution will diverge very quickly. To get around this you need to make dt smaller and/or dx and dy larger. For example, for dt=2.5e-5 and the rest as before you get alpha*dt/dx**2 + alpha*dt/dy**2 = 0.495, and the solution will look like this after 1000 iterations: Alternatively, you could use a different discretisation scheme like for ex the API scheme which is unconditionally stable but will be harder to implement.