Возникает ошибка "Индекс находился вне границ массива.", при процессе вычисления dtemp[i, j]

Код с ошибкой выделил

        int Nx, Nr, x, dt, ktp, Q;
        double R, dr, dx, a, time;
        

        Nx = 50;
        Nr = 15;
        double[,] T = new double[Nx, Nr];
        double[,] dtemp = new double[Nx, Nr];
        R = 0.15;
        x = 2;
        a = 0.00001;
        dt = 1;
        dr = R / (Nr - 1);
        dx = x / (Nx - 1);
        time = 0;

        ktp = 23; // коэффициент теплопроводности 
        Q = 2000;
        // обнуление массивов
        
        Array.Clear(T, Nx, Nr);
        
        Array.Clear(dtemp, Nx, Nr);
      

        while (time < 8000)
        {
            time = time + dt;
            time = Math.Round(time, 2);'
            foreach (int i in dtemp)
            {
                foreach (int j in dtemp)
                {
                    dtemp[i, j] = dt * a * (((T[i - 1, j] - 2 * T[i, j] + T[i + 1, j]) / Math.Pow(dx, 2)) + 1 / j * 
                        ((T[i, j] - T[i, j - 1]) / dr) + ((T[i, j - 1] - 2 * T[i, j] + T[i, j + 1]) / Math.Pow(dr, 2)));
                }
            }'
            for (int i = 0; i < Nx - 1; i++)
            {
                for (int j = 0; j < Nr - 1; j++)
                {
                    T[i, j] = T[i, j] + dtemp[i, j];
                }
            }
            for (int i = 0; i < Nx; i++)
            {
                T[i, Nr - 1] = (Q * dr) / ktp + T[i, Nr - 2];
            }

            if (time > 0)
            {
                Console.WriteLine(Math.Round(T[25,5],3));
            }
            chart1.Series[0].Points.AddXY(time, T[25, 5]);
        }

Рабочий код на Python
import numpy as np
from matplotlib import pyplot as plt

Nx=50
Nr=10

graf_x = []
graf_y = []


R=0.15
x=2

a=0.00001
dtemp=0

dt=1
dr=R/(Nr-1)
dx=x/(Nx-1)
time=0

ktp = 23 #Коэффициент теплопроводности


Q=2000

#ЗАПОЛНЯЕМ МАССИВЫ НУЛЯМИ
T=np.zeros((Nx, Nr))
dtemp=np.zeros((Nx, Nr))

#Граничные условия


#НАЧАЛО ТУТ
while time<8000: #ЦИКЛ ПО ВРЕМЕНИ
    time=time+dt
    time=round(time, 2)
    for i in range(1,Nx-1): #ЦИКЛ ПО ПРОИЗВОДНОЙ ТЕМПЕРАТУРЫ
        for j in range(1,Nr-1):
            dtemp[i,j] = dt * a * ( ((T[i-1,j]-2*T[i,j]+T[i+1,j])/dx**2) + 1/j * ((T[i,j]-T[i,j-1])/dr) \
                    + ((T[i, j-1] - 2*T[i, j] + T[i, j+1]) / dr**2) )




    for i in range(1,Nx-1): #ЦИКЛ ПО ТЕМПЕРАТУРЕ
        for j in range(1,Nr-1):
            T[i,j]=T[i,j]+dtemp[i,j]
    for i in range(0, Nx):  # ЦИКЛ ПО воздействию
        T[i,Nr-1]=(Q*dr)/ktp + T[i,Nr-2]

    #print(time)
    if time>0:
        print(round(T[25,5], 3))



    #График
    graf_x.append(time)
    graf_y.append(T[25,5])

plt.plot(graf_x, graf_y)

plt.title("Нагревание")
plt.ylabel('Temperature')
plt.xlabel('Time')
plt.show()