Testing numerical solver: Finding example Heat Equation in 2D For testing a numerical solver (FEM with linear elements with Crank Nicolson) for the heat equation with homogeneous Neumann boundary conditions.
$$\begin{cases} \frac{\partial u}{\partial t} - \Delta u = f \text{ in } \Omega \\ \frac{\partial u}{\partial n} = 0 \text{ on } \partial \Omega \end{cases}$$
on $\Omega = (0,1) \times (0,1) \subset \mathbb R^2$, where $\frac{\partial u}{\partial n}$ denotes the gradient in the direction of the normal on $\partial \Omega$, i.e. $\nabla u \cdot n$. 
I want to find suitable functions $f$ and $u$ for testing the convergence of this solver. However I wasn't able to come up with such a pair $f$ and $u$, can anyone recommend any or explain how to find examples?
 A: There's a Fourier method (separation of variables) used for homogeneous  $f = 0$ problem. Consider solution in form $u(t, x, y) = T(t) V(x, y)$
$$
T' V - T \Delta V = 0 \text{ in } \Omega\\
T \frac{\partial V}{\partial n} = 0 \text{ on }\partial \Omega
$$
Dividing the equation by $TV$ we obtain
$$
\frac{T'}{T} = \frac{\Delta V}{V}
$$
Since the left side does not depend on $x,y$ and the right does not depend on $t$, the quotient is some number $\lambda$.
Thus, $V$ should be an eigenfunction of the Laplacian with Neumann boundary conditions.
$$
\Delta V = \lambda V \text{ in } \Omega\\
\frac{\partial V}{\partial n} = 0 \text{ on }\partial \Omega\\
T(t) = T(0) e^{\lambda t}
$$
Since the problem is homogeneous, any linear combination of solutions is also a solution. I.e. you can use the following general form
$$
u(t,x,y) = \sum_{i=1}^{n} T_i(0) e^{-\lambda_i t} V_i(x, y)
$$
where $(\lambda_i, V_i)$ is an eigenpair of the Laplacian with Neumann b.c. and $T_i(0)$ are arbitrary constants. 
For the unit square domain the eigenpairs are
$$
V_{k,l}(x, y) = \cos \pi k x \cos \pi l y, \qquad k, l = 0, 1, 2, \dots\\
\lambda_{k,l} = -\pi^2 (k^2 + l^2).
$$
