Solution (source code)

= Solution

Put $\eta=\psi-\phi$. The common <Dirichlet boundary condition> gives $\eta=0$ on $S$. Expanding the <energy functional>,
$$
\begin{aligned}
E[\psi]-E[\phi]
={}&\int_V\left(|\nabla\eta|^2+m^2\eta^2\right)dV\\
&+2\int_V\left(\nabla\phi\cdot\nabla\eta
+m^2\phi\eta\right)dV.
\end{aligned}
$$
<Integration by parts> and the field equation make the cross term zero:
$$
\int_V\left(\nabla\phi\cdot\nabla\eta+m^2\phi\eta\right)dV
=\int_S\eta\frac{\partial\phi}{\partial n}\,dS
-\int_V\eta(\nabla^2\phi-m^2\phi)\,dV
=0.
$$
Consequently
$$
\boxed{
E[\psi]-E[\phi]
=\int_V\left(|\nabla\eta|^2+m^2\eta^2\right)dV\geq0}.
$$

For $m>0$, equality holds exactly when $\eta=0$. For $m=0$, equality holds exactly when $\nabla\eta=0$, so $\eta$ is constant on the connected region; its zero boundary value then again forces $\eta=0$. The minimizing property therefore proves uniqueness for the Dirichlet problem in both cases:
$$
\boxed{\psi=\phi\text{ is the equality condition for every }m\geq0}.
$$