I start here with an isolated plane and a radially symmetric heat map f. The height indicates higher temperature. I then calculate a final state function g that is the integral of f over the plane divided by the total area which is 4 in my plot.
Because the plane is isolated I figure the total thermal energy is invariant. Then I consider the homeomorphism of the heat map f transitioning to the final state map g. To simulate the fact that the closer a heat map f(t) gets to g the slower it transitions, I use a logarithmic proportionality factor A.
The final command for my simulation:
Here is the animation:
It at least looks similar to simulations of the numerical heat equation, it seems to go along with the fact that the heat flow is in the opposite direction of the temperature gradient. And it goes along well with Newton's law of cooling that each small area cools slower the closer its temperature is to it's surroundings.