Change nmax from 1 to 500.
I really like how you used meshgrid to generate all the time, position combintions as arrays Time and X, and this allows you to compute the array Tnew without using a for loop for time or a for loop for position. Nicely done!
The key was to replace matrix multiplication, "*", with element-wise multiplication, ".*" .
The plot is good now:
It is 0 degrees at x=0 and x=10, as desired. It is 100 across the bar at time t=0. It loses heat gradually with time. The temperature distribution across the bar becomes mroe rounded with time. All these results are what we expect. If you look closely at the temperature at time 0, you see "ringing" at the edges x=0 and x=10. This is the expected Gibbs phenomenon, which arises when a square wave is represented by a finite Fourier sum. If you make nmax smaller, the Gibbs phenomenon becomes more obvious.
All of this means that now you have code that computes the answer equation. But this is not a simulation. For a simulation, you should write another script, or add to this script. The new code should solve the difference equations which I derived in my previous post, for each successive time step, given the initial condition and given the boundary conditions. In this way you will fill in a surface like the surface above, starting at t=0 and progressing by with each step.
In my previous answer, I assumed that i was the index for position and j was the index for time. But your Tn is the other way around. So adjust accordingly.