from numpy import empty
import matplotlib.pyplot as pl
# Constants
L = 0.01 # Thickness of steel in meters
D = 4.25e-6 # Thermal diffusivity
N = 100 # Number of divisions in grid
a = L/N # Grid spacing
h = 1e-4 # Time-step
epsilon = h/1000
Tlo = 0.0 # Low temperature in Celsius
Tmid = 20.0 # Intermediate temperature in Celsius
Thi = 50.0 # Hi temperature in Celsius
t1 = 0.01
t2 = 0.1
t3 = 0.4
t4 = 1.0
t5 = 10.0
tend = t5 + epsilon
# Create arrays
T = empty(N+1,float) # N divisions, so N+1 points
T[0] = Thi
T[N] = Tlo
T[1:N] = Tmid
Tp = empty(N+1,float)
Tp[0] = Thi
Tp[N] = Tlo
# Main loop
t = 0.0
c = h*D/(a*a)
while t