Yes, I used fsolve. Let me ask my question more simply.
I took my pde, for which dt (really deltat) and T_e are constants,
k * del2 T = dT/dt - T_e, and rewrote it in a form,
0 = k * del2(T) - dT/dt + T_e
dT/dt is approximated using
dT/dt = (T(i) - T(i-1)) / dt, which is a kind of forward derivative.
0 = k * del2(T(i)) - (T(i) - T(i-1)) / dt + T_e
Next, I prepared a function, called heat. The variable x here is T(i) above,
the variable Tlast is T(i-1), and Te is T_e.
global K dx dy dz dt Te Tlast
z = K*del2(x,dx,dy,dz);
z = z + Te*ones(xl,yl,zl);
z = z-(x-Tlast)/dt;
y = z;
I then use this in fsolve.
Then I advance from one time to the next. It seems to work, it looks physical, and it converges well.
The question is... Is this mathematically correct? (BTW, why should one write a second order routine when Octave supplies del2?)
(My second question seems to be an adaptation of the FEM, so I will pursue that analysis as described in Wikipedia.)
Thanks for the help.