Why is the solution fleeing the graph? Logic//iteration error?

Question

Good evening. I have been writing some python code to solve the heat equation using the thomas algorithm to solve for the equations of a 2D-ADI method. The process for solving it is the following: First, i create a program which grabs a 200x320 data file (any combination of 0s and 1s will do, i use one that gives an M-shaped initial profile) and creates an initial condition for it. Therefore, i am solving over a 200x320 point mesh (x = 200, y = 320)

Then, i define the thomas algorithm, which solves the equation for all interior points( it gives me the solution U but doesen't add the Boundary conditions to it, which in this case are homogeneous BCs ( = zero). After, I program the ADI methods themselves, one which goes from the time level n to (n+1/2) and the next that goes from (n+1/2) to (n+1) , which is where we want to go from n. https://prnt.sc/wjzpp5 to see the equations I implement.

Finally, i want to represent the solution for 9 different iterations, from i = 0 to 1 = 1500.

Now, all of this is O.K ( i think) since the solution I obtain seems consistent with what i expect of a diffusion//heat equation (slowly diffusing/ "blurring"). The problem I have is that, when I graph the solution matrix, the solution just....flees the picture after some iterations!!!! (https://prnt.sc/wjzsnp, the initial profile is an M)

My suspicion here is that the problem might be in the indexes when i store the solution (thomas(a,b,c,d))after adding the boundary conditions 0 into the matrix for each iteration it goes through. Like, maybe i am not actially adding the boundary conditions correctly? or maybe i am getting my indexes wrong, so that instead of storing the thomas solution at the center of the matrix each time, i am storing the solution an index up and an index to the left ecery iteration, and as a consequence this makes it so that over time, it "leaves"?

The code i'm running is the following, any suggestion or if anyone spots where the potential problem might be are super appreciated. Apologies for my english, it is not my main language

import numpy as np
import matplotlib.pyplot as plt
from ThomasF import thomas
from matplotlib import cm

Jx = 200 #Grid w/ 200 points for x
Jy = 320 #Grid w/ 320 points for y
vx = 0.75
vy = 0.75 
#We define the thomas algorithm, which solves for the spatial points except at the boundaries
def thomas(a,b,c,d):
    n = len(b); m = n-1
    C = np.zeros(m); D = np.zeros(n)    

    C[0] = c[0]/b[0]; D[0] = d[0]/b[0]
    for i in range(1, m):
        temp = b[i] - a[i-1]*C[i-1]
        C[i] = c[i]/temp
        D[i] = (d[i] - a[i-1]*D[i-1])/temp
    D[m] = (d[m] - a[m-1]*D[m-1])/(b[m] - a[m-1]*C[m-1])

    x = np.zeros(n); x[m] = D[m]
    for i in range(1, m+1): x[m-i] = D[m-i] - C[m-i]*x[m+1-i]
    return x

#We define the function to read the initial condition file
def llegeixCI():#Reads a 200x320 datapoint file and  creates an initial condition matrix
    file_name="M_200x320.dat"
#Lectura del fitxer línia a línia
    with open(file_name) as f:
        content=f.read().splitlines()
#We define a 200x320 column matrix to store the initial condition
    un=np.zeros((200,320))
#We store the data in the matrix
    for i in range(0,200):
        j=0
    for x in content[i]:
        un[i,j]=x
        j=j+1           
return un
#To start with ( at the beginnint) our profile is the initial condition
U0 = llegeixCI()
#We define the first ADI method, which goes from n to n+1/2, and iterates over x for every y
def ADIx(vx,vy,Jx,Jy,U0):
    a = [*[vx/2] * (Jx-3)]
    b = [*[1+vx] * (Jx-2)]
    c = [*[vx/2] * (Jx-3)]
    d = zeros(Jx-2)
    U = U0
    for s in range(1,Jy-1):
        for r in range(1,Jx-1):
            d[r-1] = 0.5*vy*U[r,s+1]+(1-vy)*U[r,s]+0.5*vy*U[r,s-1]
        U[:,s] = [0,*thomas(a,b,c,d),0] #I store my results in the matrix
    return U
 #We define the second ADI method, which goes from n+1/2 to n+1, and iterates over y for every x
def ADIy(vx,vy,Jx,Jy,U0):
    a = [*[vy/2] * (Jy-3)]
    b = [*[1+vy] * (Jy-2)]
    c = [*[vy/2] * (Jy-3)]
    d = zeros(Jy-2)
    U = ADIx(vx,vy,Jx,Jy,U0)
    for r in range(1,Jx-1):
        for s in range(1,Jy-1):
            d[s-1] = 0.5*vx*U[r+1,s]+(1-vx)*U[r,s]+0.5*vx*U[r-1,s]
        U[r,:] = [0,*thomas(a,b,c,d),0]  #I store my results in the matrix
    return U 

#Graphing
i = 0
I  = iter(range(9))
prints = [0,10,50,100,200,400,700,1000,1500]


fig,ax = plt.subplots(3,3)
ax = ax.flatten()

while i<= 1500:
    if i == 0:
        U0 = llegeixCI()
    else:
         U0 =ADIy(vx,vy,Jx,Jy,U0)
    if i in prints:
        ax[next(I)].imshow(U0,cmap = cm.bone)
        #fig = plt.figure()

    i+=1