Plot MATLAB icon

Here's my attempt. To get the matrix representing the Laplacian I use LaplacianFilter on an array of symbols and CoefficientArrays to extract the coefficients.

n = 200;

shape = ArrayPad[ConstantArray[0, {n/2, n/2}], {{0, n/2}, {0, n/2}}, 1];    
shapeVector = Flatten @ Position[Flatten @ shape, 1];

symbolArray = Array[x, {n, n}];    
symbolLaplacian = LaplacianFilter[1.0 symbolArray , 1, Padding -> 0];    
lapMatrix = Last@CoefficientArrays[
    Flatten[symbolLaplacian][[shapeVector]],
    Flatten[symbolArray][[shapeVector]]];

{ev} = Eigenvectors[lapMatrix, {-1}];

result = Flatten[shape];
result[[shapeVector]] = ev;
result = Partition[result, n];

ListContourPlot[result]
ListPlot3D[result, BoxRatios -> {1, 1, 0.8}, Mesh -> False, 
  PlotStyle -> Specularity[White, 30]]

While the other answers are nice, the icon deserves a closer look:

Note, in particular, that four of the six edges are not constrained by the ostensible Dirichlet boundary conditions, nor is it clear that they solve a Neumann problem. And indeed, as I noted in the comments this is supported by the OP's first link.

In short, to produce the logo, they took 'pac-man' eigenfunctions around the convex corner, which is the trickiest bit to describe, and used a least-squares approach to impose the boundary conditions. Finally, they discarded all but two terms in the expansion:

After being so careful to satisfy the boundary conditions, the logo uses only the first two terms in the sum. This artistic license gives the edge of the logo a more interesting, curved shape.

In my view, this essentially says it is pointless to try and impose the boundary conditions at all. Since they took a certain artistic license, in replicating it one can simply try and do it by hand. We know that the plot is a function of the form $$ v(r,\theta)=\sum_{j=1,2} a_j J_\frac{2}{3}(k_j r)\sin(\tfrac{2}{3}\theta), $$ which has only four free parameters. These are enough that they can, to a good precision, simply be set by hand. I therefore used the Manipulate construct

Manipulate[
 Plot3D[(
    a BesselJ[2/3, λ Sqrt[x^2 + y^2]] Sin[
       2/3 (Arg[y - I x] + π)] + 
     b BesselJ[2/3, μ Sqrt[x^2 + y^2]] Sin[
       2/3 (Arg[y - I x] + π)]
    ) Boole[x >= 0 || y >= 0], {x, -1, 1}, {y, -1, 1},
  BoxRatios -> Automatic, Mesh -> None, Exclusions -> None, PlotPoints -> 50, Axes -> False
  ]
 , {{a, 1.275}, 0, 3}, {{b, 0.805}, 0, 3}, {{λ, 3.18}, 0, 10}, {{μ, 1.96}, 0, 10}]

to find good guesses for the parameters, which are the defaults above. The result is, I think, pretty close to the original:

One could also, if so inclined, attempt to implement their least-squares procedure. This sounds like a hazier problem to me: it is no longer "find an accurate enough solution of this specific problem", but rather "find an accurate solution that will also look like ours when you make this arbitrary restriction". It may be that many of the possible other options for the choices they make ($m$ points on the boundary, $n$ fundamental solutions) will yield similar-looking plots. But that's much too much work, I think.

I had this laying around from a course in numerical linear algebra I taught a few years ago.

Here's a matrix whose nonzero elements describe the basic shape.

size = 50;
nw = Partition[Table[i, {i, 1, size^2}], size];
sw = Partition[Table[i, {i, size^2 + 1, 2*size^2}], size];
se = Partition[Table[i, {i, 2*size^2 + 1, 3*size^2}], size];
L = ArrayFlatten[{{nw, 0}, {sw, se}}];
{m, n} = Dimensions[L];
L = Join[{Table[0, {n + 2}]},
  Map[Join[{0}, #, {0}] &, L],
  {Table[0, {n + 2}]}];

If size=4, for example, we get something like so

(* size = 4 *)
L // MatrixForm

Now, we convert this to the Laplacian.

entries[0, {_, _}] = {};
entries[k_, {i_, j_}] := Module[{goodVals},
   goodVals = Select[{
      {i + 1, j}, {i, j + 1}, {i - 1, j}, {i, j - 1}},
     L[[Sequence @@ #]] =!= 0 &];
   {k, #} -> 1 & /@ (L[[Sequence @@ #]] & /@ goodVals)];
saRules = Flatten[{Band[{1, 1}] -> -4,
    MapIndexed[entries, L, {2}]}];
lap = SparseArray[saRules];

lap is a matrix whose dimension is equal to the number of non-zero entries of L. You're essentially looking for an eigenvector of lap.

evs = Reverse[Eigenvectors[N[lap], 8,
  Method -> {"Arnoldi", "Shift" -> 0}]];

Here's a simple visualization

ev = evs[[1]];
vib = Map[If[# > 0, ev[[#]], 0] &, L, {2}];
ListPlot3D[vib, ViewPoint -> {2.3, 2.26, 1}]

OK, let's spruce it up a bit

{m, n} = Dimensions[vib];
delete[0, {i_, j_}] := If[
   (i + 1 > m || vib[[i + 1, j]] == 0 ) && 
    (j + 1 > n || vib[[i, j + 1]] == 0) &&
    (i - 1 < 1 || vib[[i - 1, j]] == 0) && 
    (j - 1 < 1 || vib[[i, j - 1]] == 0) &&
    ((i + 1 > m || j + 1 > n) || vib[[i + 1, j + 1]] == 0) &&
    ((i - 1 < 1 || j + 1 > n) || vib[[i - 1, j + 1]] == 0) &&
    ((i + 1 > m || j - 1 < 1) || vib[[i + 1, j - 1]] == 0) &&
    ((i - 1 < 1 || j - 1 < 1) || vib[[i - 1, j - 1]] == 0),
   Null, 0];
delete[x_, {_, _}] := x;
vib2 = MapIndexed[delete, vib, {2}];
ListPlot3D[vib2,
 Mesh -> None, PlotStyle -> {RGBColor[0.8, 0.2, 0.2],
   Specularity[White, 20]}, Lighting -> {
   {"Directional", RGBColor[0, 0.4, 0.4],
    {{-40, -100, 2}, {0, 0, 0}}},
   {"Directional", RGBColor[0.4, 0.4, 0],
    {{80, -50, 20}, {0, 0, 0}}}}, Background -> Black,
 Boxed -> False, ViewPoint -> {2.3, 2.26, 1}]