We have Exahausted ALL the standard Sep of Variable Shapes

Only special shapes work for separation of variables.  Any shape works with numerics.   It is surprisingly simple to compute this stuff.

Here are the evals and efuns for the unit Dirichlet Unit Disk

23.1329/4
Plot[BesselJ[0,z],{z, 0, 8}]
2.391^2

{vals,funs}=NDEigensystem[{
-Laplacian[u[x,y],{x,y}],
DirichletCondition[u[x,y]==0,True]
},
u[x,y],
 {x,y}∈Disk[{0,0}, 1/2],12];

Table[Plot3D[funs[[i]],{x,y}∈Disk[{0,0},1/2],PlotRange->All,PlotLabel->vals[[i]],PlotTheme->"Minimal"],{i,Length[vals]}]

For comparison when we computed stuff by hand we got

TableForm[
Table[ ((N[BesselJZero[n,k]])^2),{n,0, 4},{k,1,4}],
TableHeadings->{Range[0,4],Range[14]}
]

Here are the evals and efuns for the unit Dirichlet Unit Square

Method->{
"SpatialDiscretization"->{
"FiniteElement",
{"MeshOptions"->{"MaxCellMeasure"->0.0005}}
},
"VectorNormalization"->Function[{values,vectors,stiffness,damping},
s=stiffness;d=damping;
norm=vectors/Diagonal[vectors.damping.ConjugateTranspose[vectors]]^(1/2)]
}

MySquare=Rectangle[{0,0},{1,1}];
{vals,funs}=NDEigensystem[{
-Laplacian[u[x,y],{x,y}],
DirichletCondition[u[x,y]==0,True]
},
u[x,y],
 {x,y}∈MySquare,
6
];
Table[ContourPlot[funs[[i]],{x,y}∈MySquare,
PlotRange->All,
PlotLabel->vals[[i]],PlotTheme->"Minimal"],{i,Length[vals]}]

N[ 89848/(5896^2)]

d
MatrixPlot[d]

ListPlot[Eigenvalues[s,100]]

For comparison when we computed stuff by hand we got

TableForm[
Table[ (((n π))^2)+(((k π))^2),{n,1., 4},{k,1,4}],
TableHeadings->{Range[0,4],Range[14]}
]

Here are the evals and efuns for the rounded off unit Dirichlet Unit Square

p=2;
MyRoundedSquare=DiscretizeRegion[
ImplicitRegion[ ((((x-0.5))^p)+ (((y-0.5))^p))<(0.5^p), {x,y}],
Frame->True,
AspectRatio->Automatic
]

{vals,funs}=NDEigensystem[{
-Laplacian[u[x,y],{x,y}],
DirichletCondition[u[x,y]==0,True]
},
u[x,y],
 {x,y}∈MyRoundedSquare,
12
];
Table[ContourPlot[funs[[i]],{x,y}∈MyRoundedSquare,
PlotRange->All,
PlotLabel->vals[[i]],PlotTheme->"Minimal"],{i,Length[vals]}]

Now nothing can be computed by hand although we can still compare to the analytical “not rounded off” eigenvalues.

NDEigensystem