Solution (source code)

= Solution

Take a smooth <variation> $u_\varepsilon=u+\varepsilon\phi$ with $\phi=0$ on the boundary. Assume enough differentiability for the following operations, for example a smooth bounded <open set>, $u\in C^2$, and a twice continuously differentiable integrand. Differentiation under the <integral> gives
$$
\left.\frac{dI[u_\varepsilon]}{d\varepsilon}\right|_{0}=\int_{\mathcal D}(L_u\phi+L_{u_x}\phi_x+L_{u_y}\phi_y)\,dx\,dy.
$$
By <integration by parts>, the boundary term vanishes and the <first variation> is
$$
\int_{\mathcal D}\left(L_u-\partial_xL_{u_x}-\partial_yL_{u_y}\right)\phi\,dx\,dy.
$$
Stationarity for every such <variation>, including every compactly supported smooth one, and the <fundamental lemma of the calculus of variations> give the <Euler-Lagrange equation>
$$
\boxed{L_u-\partial_xL_{u_x}-\partial_yL_{u_y}=0}.
$$
For the quadratic integrand, $L_u=2k^2u$, $L_{u_x}=2u_x$, and $L_{u_y}=2u_y$. Thus
$$
\boxed{u_{xx}+u_{yy}=k^2u}
$$
with the prescribed <Dirichlet boundary condition>. For real $k$, a stationary function is also the unique minimizer, if it exists: for any nonzero admissible $\phi$, the linear term vanishes and $I[u+\phi]-I[u]=\int_{\mathcal D}(|\nabla\phi|^2+k^2\phi^2)>0$. This last positivity uses the zero boundary values; it also holds at $k=0$.