¿Cómo puedo optimizar un código que utiliza arrays de forma intensiva en Python?

Disculpa, esta es probablemente una pregunta muy básica, pero estoy convirtiendo algún código que he estado modelando desde MATLAB a Python, tanto para ayudarme a aprender Python como para ver si podría ejecutarlo más rápido. En MATLAB, este código tarda aproximadamente 1 segundo en ejecutarse, pero en Python, tarda alrededor de 1 minuto. ¿Hay alguna forma de acelerarlo, o esta no es una buena aplicación de Python?

import numpy as np
import matplotlib.pyplot as plt

N = 7e5 #cantidad de pasos de tiempo
dt = 1e-6 #Paso de tiempo, en segundos
tf = dt*N #Tiempo final, segundos.
trange = np.linspace(0,tf,int(N+1)) #rango de tiempo

dx = L/M #tamaño de paso espacial en termoeléctrica, metros

<h1>Define número de Fourier adimensional en termoeléctrico</h1>

Fo = dt(k/c)/(dx*2)

<h1>perfil de temperatura en termoeléctrica como función de espacio y tiempo</h1>

T = np.zeros((M+1,2))

<h1>Asignar condición inicial</h1>

T[:,0] = Ti

<h1>Establecer condición de frontera en x=L</h1>

T[M,:] = T2

<h1>perfil de temperatura frente al tiempo del lado frío de la termoeléctrica</h1>

coldTemp = np.zeros(len(trange))

<h1>temperatura inicial en lado frío</h1>

coldTemp[0] = Ti

<h1>ajustar la corriente al valor óptimo de corriente continua</h1>

I = Issmax

<h1>iterar sobre pasos de tiempo</h1>

for p in range(int(N)):

#Usar el método central de diferencias en el tiempo para encontrar la temperatura dentro
#del material termoeléctrico.
for n in range(M-1):
    #calcular cambio de temperatura en el siguiente paso de tiempo
    T[n+1,1] = T[n+1,0] + Fo*(T[n+2,0]-2*T[n+1,0]+T[n,0]) + dt*((I)**2*rho/(c*d**2*w**2))

#Aplicar el balance de energía al metal (supuesto isotérmico) y usar el
#hecho de que la temperatura del metal es igual a la temperatura del termoeléctrico
T[0,1] = T[0,0] + dt*((I)**2*rhom/(cm*lm**2*wm**2)) - (k*dt/(cm*dx*lm))*(T[0,0]-T[1,0]) - (dt*(I)*S*T[0,0]/(cm*d*w*lm))
#Guardar temperatura en el lado frío
coldTemp[p+1] = T[0,1]
#Establecer el perfil de temperatura actualizado
T[:,0] = T[:,1]

Graficar temperatura en el lado frío vs tiempo

plt.plot(trange, coldTemp)