Define the two-dimensional fundamental solution:

g[x_,y_,t_,k_]=Exp[-(x^2+y^2)/4/k/t]/4/Pi/k/t

A first plot:

Plot3D[g[x,y,1,1],{x,-5,5},{y,-5,5},PlotRange->All]

Improve the plot by using a z=constant mesh rather than x=constant, y=constant mesh and by increasing the number of points at which the function is evaluated:

Plot3D[g[x,y,1,1],{x,-5,5},{y,-5,5},PlotRange->All,
MeshFunctions->(#3&),
PlotPoints->50]

Animate:

Animate[
Plot3D[g[x,y,t,1],{x,-5,5},{y,-5,5},PlotRange->{0,1},
MeshFunctions->(#3&),
PlotPoints->25],
{t,0.001,1},
AnimationRunning->False]

The animation is not great because the time to compute each frame is long so Mathematica does a rough version to keep up.  In order to get a better animation, we can generate all of the frames in advance and then animate.  To generate the frames, we use the Table command:

frames=Table[
Plot3D[g[x,y,t,1],{x,-5,5},{y,-5,5},PlotRange->{0,1},
MeshFunctions->(#3&),
PlotPoints->50],
{t,0.001,1,0.01}];

We ended the command above with a semi-colon (;).  This suppresses the display of output since we don’t need to see all 100 frames.  Now we can use ListAnimate to play the frames in an animation:

ListAnimate[frames,AnimationRunning->False]