파이썬으로 이송확산 방정식을 통한 수치해석에 대해 공부하고 있는 학생입니다.. 자꾸 코드를 실행하면 c와 rhs배열에 nan이 들어가는데 그 이유를 모르겠습니다..

#Library import
import numpy as np
import math
import matplotlib.pyplot as plt

#Definitions and initializations of variable arrays
m = 100
v = 1.0
dx = 1.0
dt = 0.02
nmax = 2250 # T / dt (T=45)
x = np.zeros(shape = m+1) # 0~100 array 
c = np.zeros(shape = m+1)
rhs = np.zeros(shape = m+1) # => 행벡터
a1 = np.zeros(shape = m+1)
b1 = np.zeros(shape = m+1)
c1 = np.zeros(shape = m+1)
d = 0.0
e = 0.0
f = 0.0
c0 = np.zeros(shape = m+1)
ce = np.zeros(shape = m+1)

# Selection of scheme
alpha = float(input('input alpha : ')) #키보드 입력 , cf) 강제 형변환 
theta = float (input('input theta : ')) #0, 0.5, 1.0  (theta가 1일때 implicit method, 0일때 explicit method)
pe = float(input('input pe : '))#0.1 , 1
scheme = float(input('input scheme :')) 

#Grid generation
#% when i = 0, its distance is at (1*dx)m from origin
for i in range (0, m+1):
    x[i] = i*dx
    
#initial condition
for i in range (0, m+1):
    if i == 0:
        c[i] = 1.0
    else :
        c[i] = 0.0

#Save the initial condition
for i in range (0, m+1):
    c0[i]=c[i]
       
#==========================================
#Main Program
#==========================================
#Implicit method 
#Time integral

for iter in range (1,nmax+1): #for 1문
    if iter == round(iter/100)*100: #100번 iteration마다 찍고 싶을 때 사용
        print('{0}'.format(iter))
    for i in range (1,m): ## for 2문
        a1[i] = (-v*alpha*theta/(2*dx))-v*(1-alpha)*theta/dx-v*theta/(pe*dx);
        b1[i] = 1/dt + (v*theta*(1-alpha)/dx) + (2*v*theta/(pe*dx))
        c1[i] = (v*alpha*theta/(2*dx))-(v*theta/(pe*dx))
        #RHS Terms 
        d = (v*alpha*(1-theta)/(2*dx))+(v*(1-alpha)*(1-theta)/dx)+(v*(1-theta)/(pe*dx))
        e = 1/dt - (v*(1-alpha)*(1-theta)/dx)-(2*v*(1-theta)/(pe*dx))
        f = (-v*alpha*(1-theta)/(2*dx)) + (v*(1-theta)/(pe*dx))
        rhs[i] = d*c[i-1]+e*c[i]+f*c[i+1] ##end of for 2문
        
    rhs[1]=rhs[1]-a1[1]*c[0]
    rhs[m-1]=rhs[m-1]-c1[m-1]*c[m]
         
#Thomas algorithm
    for i in range(2,m): ##for 3문
        r=a1[i]/b1[i-1]
        b1[i]=b1[i]-r*c1[i-1]
        rhs[i]=rhs[i]-r*rhs[i-1] ##end of for 3문
     
    rhs[m-1]=rhs[m-1]/b1[m-1];
    for i in range(m-2,0,-1): ##for 4문
        rhs[i]=(rhs[i]-c1[i]*rhs[i+1]/b1[i])##end of for 4문
    #end of Thomas algorithm
    
    #renewal of concentration at the next time step
    for i in range (1, m): ##for 5문
        c[i]=rhs[i] ## end of for 5문
    #renwal of boundary conditions 
    c[0] = 1.0 #continous source injection
    c[m] = 0.0 #Neumann boundary condition
    # end of for 1문
    
#Exact Solution
for i in range (0, m+1):
    d = v*dx/pe
    y1 = (x[i]-v*dt*nmax)/(2*math.sqrt(d*dt*nmax))
    y2 = (x[i]+v*dt*nmax)/(2*math.sqrt(d*dt*nmax))
    ce[i]=c[0]/2*(math.erfc(y1)+math.exp(v*x[i]/d)*math.erfc(y2))