Проблема с циклом Python: итерация конечного объема/разности

Мне трудно понять, почему этот цикл for не работает. Полный код также приведен ниже.

Я хочу рассчитать температуру в будущие моменты времени (j+1) на основе температур в предыдущий момент времени (j), используя метод конечных объемов. У меня есть функции настройки для расчета входной и выходной проводимости, а также внутренней генерации внутри узла. У меня есть еще одна функция, которая суммирует их и определяет результирующую температуру на основе уравнения энергетического баланса. Однако моя проблема связана именно с циклом ниже:

for j in range(1, n): # time
    for e in range(1, i-1): # position
        T[e,j] = Tn_plus_1( T[e, j-1] , timestep,rho[e], cp[e], re_nodes[e] , rw_nodes[e],
                
                conduction_in(k[e-1], rw_nodes[e],T[e, j-1], T[e-1, j-1], r_nodes[e], r_nodes[e-1]),
                
                conduction_out(k[e], re_nodes[e], T[e+1, j-1], T[e, j-1], r_nodes[e+1], r_nodes[e] ),
                
                fraction[j]*internal_gen(q_dens[e], timestep, rho[e], cp[e]) )

Tn_plus_1 — это функция для расчета температуры узла на следующем временном шаге на основе энергетического баланса и метода конечного объема. Входная, внешняя проводимость и внутренняя генерация — это функции для расчета теплопередачи каждого механизма в конкретный узел или из него. Я настроил циклы для итерации по времени и положению, чтобы определить изменение температуры в каждом узле на протяжении всего анализа.

Проблема в том, что вывод проводимости неправильно применяется к узлу 7. Я подозреваю, что та же проблема возникает с выводом проводимости во втором узле. Мыслительный процесс, лежащий в основе цикла, заключается в том, что я хочу вычислить температуры внутренних узлов и сохранить граничные условия на двух внешних узлах. Следовательно, я устанавливаю цикл между 1 и i-1 (где 0 и i являются границами).

Всего существует 9 узлов положения (каждый из r, re и rw имеет по 9 значений). Также имеется 9 значений свойств материала (rho, cp, k).

Я новичок в Python и подозреваю, что проблема в том, что я запутался с индексацией Python, но я не понимаю, где именно цикл идет не так для двух узлов рядом с каждой границей. Заранее спасибо.

import numpy as np
import matplotlib.pylab as plt

empty = np.zeros(9) # empty array to fill with required properties

#material properties (will vary with position at a later stage)
k = empty + 23
rho = empty + 1600
cp = empty + 1700

#node spacing
dr_internal = 0.15
dr_external = 0.14

r_nodes = np.linspace(0, 1.2, 9) # nodes, m

r_nodes2 = np.insert(r_nodes,9 , r_nodes[-1]+dr_external) #inserting additional point to determine
east boundary on last node
re_nodes = (r_nodes2[1:] + r_nodes2[:-1]) / 2 #east locations
rw_nodes = np.insert(re_nodes, 0, 0)[:-1] #west locations

#q density for the internal generation term
q_dens = empty + 4400000
q_dens[0] = 0
q_dens[1] = q_dens[1]*0.625
q_dens[-1] = q_dens[-1]*0.5

#setting up time increments for the analysis
timestep = 15
steps = 100
total_time = timestep * steps
time = np.arange(0, steps*timestep, timestep)

#defining the number of time points and spacial nodes
i = len(r_nodes)
n = len(time)
T = np.zeros((i,n))

#initialising temperatures based on the boundary condition and initial condition
IC = 873
BC = [873, 873]
T[0,:] = BC[0]
T[-1,:] = BC[1]
T[:,0] = IC

#energy generation is time dependent, so creating a fraction to multiply internal generation term by
fraction=np.zeros(steps)
for j in range(1, steps):
    time2 = j*timestep
    fraction[j] = 0.06*time2**(-0.2)
fraction[0] = 0.06


def conduction_in( kw, rw, Ti, Ti_minus_1, ri, ri_minus_1):
    return (-kw*rw*((Ti-Ti_minus_1)/(ri-ri_minus_1)))

def conduction_out(k, re, Ti_plus_1, Ti, ri_plus_1, ri):
    return (k*re*((Ti_plus_1-Ti)/(ri_plus_1-ri)))

def internal_gen(q_dens, dt, rho, cp):
    return q_dens*dt/(rho*cp)

def Tn_plus_1(Ti, dt, rho, cp, re, rw, conduction_in, conduction_out, internal_gen):    
        return Ti + ((2*dt/(rho*cp*(re**2-rw**2)))*(conduction_in - conduction_out ))+ internal_gen

for j in range(1, n): # time
    for e in range(1, i-1): # position
        T[e,j] = Tn_plus_1( T[e, j-1] , timestep,rho[e], cp[e], re_nodes[e] , rw_nodes[e],
                
                conduction_in(k[e-1], rw_nodes[e],T[e, j-1], T[e-1, j-1], r_nodes[e], r_nodes[e-1]),
                
                conduction_out(k[e], re_nodes[e], T[e+1, j-1], T[e, j-1], r_nodes[e+1], r_nodes[e] ),
                
                fraction[j]*internal_gen(q_dens[e], timestep, rho[e], cp[e]) )
        

plt.plot(T[:,99]) 

Что я получаю:

Чего я ожидаю: