T_num = nan(length(z),length(t));
T_num(:,1) = T_av;
T_num(1,n+1)=-sur_flux*dz/conductivity + 250;
for j = 2:J-1
T_num(J,n+1) = -flux*dz/conductivity+T_num(J-1,n+1);
The above first sets first column to T_av
The loop on sur_flux isn't what you intend I think; it sets the same location of (1,n+1) sequentially to the computed result but does that without changing n so all that happens before the j loop is to set
T_num(1,n+1)=-1350*dz/conductivity + 250;
all the other values are immaterial because the same node is overwritten every pass. Not exactly sure what your intent is, but that undoubtedly isn't it... :)
Then with n=1, for j=2,
Hmmmm...that looks ok; I thought there was going to be a reference to an uninitialized location that was going to propagate through.
I still think that must be so and I missed it...set a breakpoint and step through with the debugger and see where the NaN comes from...