I have a problem with coding a Diffusion equation (FT-CS method)

Question

In statistical physics, we encounter problems that we can express via differential equations, first order, second order, time based, space based .... and there are some problems that are expressed throughout both time and space, in this case, heat diffusion.

what I am trying to do in to simulate a heat transition from a hot heat source to a cold heat source via a thermal conductor, in this case a metal bowl at 20°C, immersed in a cold bath 0°C, filled with hot water 50°C.

this is the code (I tried Chat GPT to figure out the problem) :

import numpy as np
import matplotlib.pyplot as plt

# Parameters
m = 50  # Number of spatial grid points
D = 4.25e-6  # Diffusion coefficient
L = 1.0  # Length of the domain
T = 1.0  # Total simulation time
delta_x = L / m  # Spatial step size
delta_t = 0.4 * delta_x**2 / D  # Time step size (adjusted for stability)

# Initialize temperature array
phi = np.ones((m+2, m+2)) * 20  # Initialize with a uniform temperature of 20°C
phi[:, 0] = 50  # Boundary condition at y = 0 (left boundary)
phi[:, -1] = 0  # Boundary condition at y = L (right boundary)

# Plot initial temperature distribution as a map
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.imshow(phi, origin='lower', extent=[0, L, 0, L], cmap='jet', aspect='auto')
plt.colorbar(label='Temperature (°C)')
plt.xlabel('Position (x)')
plt.ylabel('Position (y)')
plt.title('Initial Temperature Distribution')

# Iteration
num_steps = int(T / delta_t)
for step in range(num_steps):
    phi_old = np.copy(phi)  # Store the current temperature distribution
    for i in range(1, m+1):
        for j in range(1, m+1):
            # Update temperature using the FTCS method
            phi[i, j] = phi_old[i, j] + D * delta_t / delta_x**2 * (
                phi_old[i+1, j] + phi_old[i-1, j] + phi_old[i, j+1] + phi_old[i, j-1] - 4*phi_old[i, j])
            
            # Apply boundary conditions
            phi[i, 0] = 50  # Boundary condition at y = 0 (left boundary)
            phi[i, -1] = 0  # Boundary condition at y = L (right boundary)

# Plot final temperature distribution as a map
plt.subplot(1, 3, 2)
plt.imshow(phi, origin='lower', extent=[0, L, 0, L], cmap='jet', aspect='auto')
plt.colorbar(label='Temperature (°C)')
plt.xlabel('Position (x)')
plt.ylabel('Position (y)')
plt.title('Final Temperature Distribution')
# Plot final temperature distribution as a curve
plt.subplot(1, 3, 3)
plt.plot(np.linspace(0, L, m+2), phi[:, int((m+1)/2)], color='blue')
plt.xlabel('Position (x)')
plt.ylabel('Temperature (°C)')
plt.title('Temperature Distribution along x-axis')

plt.tight_layout()
plt.show()enter image description here

the final code is supposed to show a heat map of some sort for the initial condition, final result, and a curve of temperature and position.

the final result is supposed to be a gradient from hot to cold, and a linear-ish curve.

the only problem is that the final results are identical to the initial result, even though physically it was supposed to be equilibrium at 25°C, it is still 20, and the curve is still 50°C top, 0°C bottom, and 20°C all around.

I did everything, tried to plot with phi_old instead of phi, checked the steps, checked the boundary enter image description here conditions are met inside the loop, nothing works.