I am just learning PDE, and I am interested to compare analytical solution with numerical solution of Helmholtz equation in a unit square with zero boundary condition. I am not sure if it possible. Below is the equation I am analyzing, with its eigenvalues and eigenfunctions:
$$
\begin{cases}
\Delta^2u+\lambda u=0,\\
u(0,y)=0,\quad u(x,0)=0,\\
u(1,y)=0,\quad u(x,1)=0,
\end{cases}
$$
I found out from a textbook that the eigenvalues and eigenfunctions are:
$\lambda_{nm}=(n^2+m^2)\pi^2$ and $u_{nm}=\sin n\pi x \sin m\pi y$, where $n,m=1,2,3,\ldots$.
I know that the solution to the above equation is trivial, i.e., $u=0$. However, I found from question regarding solution to Helmholtz by other user, that we can solve the equation in term of eigenfunctions and eigenvalues. So, in my case, the first six eigenfunctions are 
region = Rectangle[];
{eigenvalues[region], eigenfunctions[region]} = NDEigensystem[{-Laplacian[u[x, y], {x, y}] + u[x, y], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} \[Element] region, 6];
Grid[Partition[Table[Show[{ContourPlot[eigenfunctions[region][[j]], {x, y} \[Element] region, Frame -> None, ColorFunction -> "BlueGreenYellow", PlotPoints -> 60, PlotRange -> Full, PlotLabel -> eigenvalues[region][[j]]],RegionPlot[region, PlotStyle -> None, BoundaryStyle -> {Black, Thick}]}], {j, 1, Length[eigenvalues[region]]}], 3]]
.
I wonder if I can do error analysis by suming the eigenfunctions and compare it with $u=0$ in $L_2$ norm.
But, what is the exact analytical solution of this equation in terms of eigenfunctions? I know it is the summation of those eigenvalues, but I am not sure how to sum it.
And, how to sum those first six images of eigenfunctions with mathematica? Once somebody show me how to sum those images, I think I can do the error analysis myself. Thanks for your time!