FDTD ecuacion de difusion de calor

import numpy as np
import matplotlib.pyplot as plt
import math
import cmath


kappa_al  = 237.0;
rho_al    = 2698.4;
c_al      = 900.0;
D         = kappa_al/(rho_al*c_al);

Nx        = 10
xi        = -0.04;
xf        =  0.04;
dx        = (xf-xi)/Nx;
x         = np.zeros(Nx+1)
TTa       = np.zeros(Nx+1)
TTn       = np.zeros(Nx+1)
t0        = 0.1

Nt        = 100;
dt        = 0.2*(dx*dx)/(2.0*D);
t         = np.zeros(Nt+1)

#
it        = 0
t[it]     = t0

for ix in range(0,Nx+1):
  x[ix]  = xi+ix*dx
  TTa[ix]= math.sqrt(t0/t[it])*math.exp(-(x[ix]*x[ix])/(4.0*D*t[it]))
  TTn[ix]= TTa[ix]
    
plt.ylim(0,1)
plt.plot(x,TTa,'-ob')
plt.savefig("T_00{}.png".format(0))
plt.close()

for it in range(1,Nt+1):
  t[it]     = t0+it*dt
  for ix in range(1,Nx):
    TTa[ix]= math.sqrt(t0/t[it])*math.exp(-(x[ix]*x[ix])/(4.0*D*t[it]))
    TTn[ix]= TTn[ix] + ((D*dt)/(dx*dx))*(TTn[ix+1]-2.0*TTn[ix]+TTn[ix-1])

  plt.ylim(0,1)
  plt.plot(x,TTa,'-ob',x,TTn,'-or')
  if it < 10:              plt.savefig("T_00{}.png".format(it))
  if it >= 10  and it<100: plt.savefig( "T_0{}.png".format(it))
  if it >= 100 and it<100: plt.savefig(  "T_{}.png".format(it))
  plt.close()