Since we are not four dimensional (or more) beings, it is hard four us to look inside three dimensional objects. Unless, of course, they have some kind of openings. An easy way to open up a cube is to make it into lots of small cubes, like this.
It is just a lot of Cuboid's. The color function which I called heatColor is essentially a wrapper around RGBColor. It has two parameters min and max, for temperatures T close to min we choose litte red and a lot of blue. Vice versa for T close to max.
In order to solve the partial differential equation I use NDSolve. The boundary conditions are zero on all sides of the cube except the bottom (z=0). I just chose something that would be consistent with the other boundary conditions as well as with the inital condition. The inital condition is also chosen as anything simple that will be zero when x, y or z get close to zero or one. The result from NDSolve comes in the form of an interpolated function, so you can evaluate the numerical solution in any point that you wish!
You can play around with the number of cubes and the opacity to find your personal optimum.
The code:
solution =
NDSolve[{D[u[x, y, z, t], t] ==
D[u[x, y, z, t], x, x] + D[u[x, y, z, t], y, y] +
D[u[x, y, z, t], z, z],
u[0, y, z, t] == 0,
u[1, y, z, t] == 0,
u[x, 0, z, t] == 0,
u[x, 1, z, t] == 0,
u[x, y, 0, t] == t x y (1 - x) (1 - y),
u[x, y, 1, t] == 0,
u[x, y, z, 0] == x y z (1 - x) (1 - y) (1 - z)},
u,
{x, 0, 1}, {y, 0, 1}, {z, 0, 1}, {t, 0, 2}]
(* Let v be the first - and only - solution *)
v = (u /. solution)[[1]]
heatColor[T_, min_, max_] :=
With[{c = (T - min)/(max - min)}, RGBColor[{c, 0, 1 - c}]]
Graphics3D[Table[
{heatColor[v[i/5, j/5, k/5, 1], 0, 0.03],
Opacity[.50],
Cuboid[{i, j, k}/5, {i, j, k}/5 + .10]},
{i, 0, 5}, {j, 0, 5}, {k, 0, 5}], Axes -> True,
AxesLabel -> {x, y, z}]
