I numerically solved the heat equation using an implicit scheme. The solution converges, however, it bothers me that the values of the function from layer to layer of the grid change sign. Is such behavior possible or did I make a mistake?
I attach the part of the values of the function in the nodes of the grid.
Certainly having negative heat content is suspicious.
My first thought (which turned out to be wrong) was that your time step size was too large. In places where the second derivative in space is negative, the downward update applied for too long can result in negative values. This is not what is going on in your computation.
It really just looks like you have (or need) an extra minus sign on the output of your update rule.