import numpy as np
import matplotlib.pyplot as plt

Lambda = 0.037
Cp = 1500
Rho = 1.325
L = 1
t_max = 20000
N_t = 100
N_X = 5
T_int = 20
T_ext = 5
K = Lambda/(Rho*Cp)

dx = L/N_X
dt = t_max/N_t
Temp = np.zeros((N_t+1,N_X+1))

Temp[0,0] = T_int

for i in range (1,N_X+1):
    Temp[0,i] = T_ext
    
for n in range (1,N_t+1):
    Temp[n,0] = T_int
    Temp[n,N_X] = T_ext
    
for n in range (0,N_t):
    for i in range(1,N_X):
        Temp[n+1,i] = (dt*K/dx**2)*(Temp[n,i+1]+Temp[n,i-1]-2*Temp[n,i]) + Temp[n,i]

x = np.linspace(0,L,N_X+1)  

plt.plot(x,Temp[0,:],'-o')
plt.plot(x,Temp[33,:],'-')
plt.plot(x,Temp[66,:],'-+')
plt.plot(x,Temp[100,:],'--')

plt.show()