Weak form of the PDE $-\nabla^2u -k^2 u=0$ I'm new to this weak formulation theories, and I'm having a hard time deriving the weak form to following PDE:
$$E(u)=-\nabla^2u -k^2 u=0$$
using the fact that for the minimizer $\displaystyle\lim_{\epsilon\to 0}\frac{d}{d\epsilon}E(u+\epsilon v) = 0$


*

*Is this method general for all kind of PDE equations (integral form and ...)? see link below
Derivation of weak form for variational problem

*If not, could you suggest any other method? Thanks.

 A: 
Is this method general for all kind of PDE equations (integral form and ...)? see link below:Derivation of weak form for variational problem

This method actually originates from calculus of variations when looking for the Euler-Lagrange equation for a functional. This functional usually describes certain potential and kinetic energy (this is why people usually call this functional energy functional). For example, from electrostatic potential energy functional to Poisson equation. The expression 
$$
\lim_{\epsilon\to 0} \frac{d}{d\epsilon}E(u+\epsilon v) = 0 \;\text{ for any }v\in V \iff u = \operatorname{argmin}\limits_{v\in H} E(v).
$$


$E(u)=-\nabla^2u -k^2 u=0.$

This is not accurate. $E(u)$ should be identified with the energy functional I mentioned above, not the underlying PDE you gave, which is the Helmholtz equation (or say the eigenvalue problem for $-\Delta = -\nabla^2$)
$$
-\Delta u = k^2u.\tag{1}
$$
The energy functional you wanna minimize is actually the Rayleigh quotient:
$$
E(v) = \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2},
$$
given the Rayleigh quotient is $k^2$. The weak form you will obtain will be (left for you as an exercise, this is called the Minimum Principle of the Eigenvalue)
$$
\lim_{\epsilon\to 0} \frac{d}{d\epsilon} E(u+\epsilon v) = 0 \implies \int_{\Omega} \nabla u\cdot \nabla v - k^2 \int_{\Omega} uv = 0.\tag{2}
$$


If not, could you suggest any other method? Thanks.

A simpler, but containing less physics based heuristics, method is that multiplying both sides of (1) by a test function $v$, then integration by parts using Green's formula (or divergence theorem, Stokes, etc):
$$
\int_{\Omega }-\Delta u v =k^2\int_{\Omega } uv \implies  \int_{\Omega} \nabla u\cdot \nabla v -\color{blue}{\int_{\partial \Omega} v\frac{\partial u}{\partial \nu} \,dS}= k^2 \int_{\Omega} uv.
$$
Then we set the test function $v=0$ on the boundary to make blue term vanish. Hence you have (2). The reason to do this for Dirichlet eigenvalue problem (equation (1) with Dirichlet boundary value) is that: The test functions should not change the boundary value of the solution to the underlying PDE.

Some remark: in the question you link, the OP's notation was a bit inaccurate, the bilinear form should NOT be identified with the same notation with the energy functional.
