Illustrating Solutions to 2–D Partial Differential Equations

Mathematica code as it appeared in Mathematica in Education 
Volume 1 Number 1, Fall 1991.

by Matthew M. Thomas
Department of Chemical Engineering
Washington University in St. Louis
thomas@wuche1.wustl.edu

Figures 1 & 2

W[x_, y_, t_] := (64/N[Pi]^6) *
                 N[Sum[(Cos[N[Pi] t Sqrt[m^2 + n^2]] -
                       (Sin[N[Pi] t Sqrt[m^2 + n^2]]/
                       (N[Pi] Sqrt[m^2 + n^2])))*
                 Sin[m N[Pi] x] Sin[n N[Pi] y]/(m^3 n^3),
                 {m, 1, 15, 2}, {n, 1, 15, 2}]];

GrayPlot[f_, ranges__, n_] :=
   Plot3D[f, ranges,
      Mesh -> True,
      PlotRange -> {-0.07, 0.07},
      PlotPoints -> n];

WavePlot[n_, t_] :=
   GrayPlot[W[x, y, t], {x, 0, 1}, {y, 0, 1}, n]

WavePlot[20, 0.0]

Figure 3

U[x_, t_] := Block[{k=Ceiling[1/Sqrt[t]]},
             2 Sum[Block[{nu=(m+.5) N[Pi]},
             (-1)^m/nu Exp[-nu^2 t] Cos[nu x]],
                  {m, 0, k}]];

U[x_, 0.] := 1;

ListHuePlot[array_] :=
    ListPlot3D[
        Table[Table[0, {Length[array]+1}], {Length[array]+1}],
        Map[GrayLevel[.75 #]&, array, {2}],
        Boxed -> False, ViewPoint -> {0, 0, 10},
        Mesh -> False, Axes -> None];

HeatPlot[n_, t_] := Block[{u=Table[U[x, t], {x, -1, 1, 2/n}]},
    ListHuePlot[Outer[Times, u, u]]];

HeatPlot[50,0.1]   (* to produce Fig. 3 densityplot *)