For a start I think you’ve missed declaring dx2 somewhere in your code. Check out my comment below:
import numpy as np
import matplotlib.pyplot as plt
L = 10
n = 10
T0= 0
T_rock = 600
T_magma = 1300
alpha = 10**(-6)
dx = L/n
t_final = 60
dt = 0.1
x = np.linspace(dx/2, L-dx/2, n)
t = np.arange(0, t_final, dt)
dTdt = np.empty(n)
T = np.ones(n)*T0
for j in range(1, len(t)):
for i in range(1, n-1):
dTdt[i] = alpha*(-(T[i]-T[i-1])/dx2 + (T[i+1]-T[i])/dx2) # <---- where did dx2 come from?
dTdt[0] = alpha*(-(T0-T_magma)/dx2 + (T[1]-T0)/dx2)
dTdt[n-1] = alpha*(-(T[n-1]-T[n-2])/dx2 + (T_rock-T[n-1])/dx2)
T = T + dTdt*dt
plt.figure(1)
plt.plot(x,T)
plt.axis([0, L , 500, 1500])
plt.xlabel('Distance [m]')
plt.ylabel('Temperature [°C]')
plt.pause(0.05)
plt.show()
Please post your traceback in tripple back tics for any subsequent errors.
They’re really helpful for us trying to help you 