Laplace equation from a functional solved on a disk - arising from a functional equation (Note this is a completely updated version of the original question with a different approach).
I have the given problem:

Let $\Omega$ be the disk $\Omega:x^2+y^2<R^2\\$
Determine all extremals of the functional:
$J[z]=\int\int_\Omega\big((z_x)^2+(z_y^2)\big)dxdy, \ \ \ z=z(x,y)$,
under the boundary condition that $z=0$ along the curve
$\partial\Omega$ and the additional  condition that
$\int\int_\Omega z^2dxdy=1$

Since this  functional yields the Laplace equation from using the following formula
\begin{equation}
F_z-\frac{\partial}{\partial x}F_{z_x}-\frac{\partial}{\partial y}F_{z_y}=0
\end{equation}
we have to consider the problem as a  Laplace equation on radial coordinates.
So the problem to solve is:
\begin{equation}
z_{xx}+z_{yy}=\lambda z , \ \ \ \Omega:x^2+y^2\leq R^2, \\
\partial\Omega=0
\end{equation}
So these are Dirichlet conditions on a disk, which we write as a rectangle in polar coordinates: $\{x^2+y^2\leq R^2\}=[0,R)\times[0,2\pi)$ and shall solve with the variable change $z(x,y)\rightarrow z(r,\phi)$
We therefore solve the following Dirichlet problem in polar coordinates by separation of variables
\begin{equation}
z_{rr}+\frac{1}{r}z_r+\frac{1}{r^2}z_{\phi\phi}=0
\end{equation}
By separation of variables, we obtain two ODEs:
\begin{equation}
R_{rr}+\frac{1}{r}R_r-\lambda R=0\\
\Phi_{\phi\phi}=-\lambda\Phi
\end{equation}
Solving each separately, we obtain
\begin{equation}
R(r)=ar^\lambda+br^\lambda\\
\Phi(\phi)=A_n\cos \lambda\phi+B_n\sin \lambda\phi
\end{equation}
The eigenvalue $\lambda=n$ and Dirichlet conditions give $b=A_n=0$:
\begin{equation}
R(r)=ar^n\\
\Phi(\phi)=B_n\sin n\phi
\end{equation}
So we have the general solution, where $\beta_n=a\cdot B_n$:
\begin{equation}
z(r\phi)=\frac{A_0}{2}+\sum_{n=1}^\infty r^n\beta_n\sin n\phi
\end{equation}
Now, since the rectangle converted to a disk we re-write the following Fourier coefficients:
\begin{equation}
A_0=\frac{1}{L}\int_0^L f(x)dx, \text{and}\\
\beta_n=\frac{2}{L}\int_0^Lf(x)\sin\frac{k\pi x}{L}dx
\end{equation}
as:
\begin{equation}
A_0=\frac{1}{\sqrt{\pi R^n}}\int_0^{2\pi} f(\phi)d\phi, \text{and}\\
\beta_n=\frac{2}{\sqrt{\pi R^n}}\int_0^{2\pi}f(\phi)\sin n\phi d\phi
\end{equation}
We now have the full form of the solution:
\begin{equation}
z(r\phi)=\frac{1}{2\sqrt{\pi R^n}}\int_0^{2\pi} f(\phi)d\phi+\sum_{n=1}^\infty \frac{2}{\sqrt{\pi R^n}}\int_0^{2\pi}f(x)\sin n\phi d\phi\sin n\phi
\end{equation}
However, we have been given Dirichlet conditions, without any initial conditions, $z(x,0)=f(x)$ . But instead, we have been given the IC:
\begin{equation}
\int\int_\Omega z^2dxdy=1
\end{equation}
which is rather useful. So we write up the integral in polar coordinates with the necessary information for $z$, where the coefficient $\frac{A_0}{2}$ is not relevant, since the IC.s are given by that simple double integral which we rewrite to polar form:
\begin{equation}
\int\int_\Omega z^2dxdy=1=\int\int_\Omega z^2rdrd\phi=1
\end{equation}
where $\Omega_{x,y}\rightarrow\Omega_{r\phi}$
This gives simply
\begin{equation}
\int_0^{2\pi}\int_0^R\bigg[\sum_{n=1}^\infty r^n\beta_n\sin n\phi\bigg]^2rdrd\phi=1
\end{equation}
when I solved this, I got:
\begin{equation}
\beta_n=\sqrt{\frac{(4n+2)}{2\pi R^{2n+1}}}
\end{equation}
So the final solution is:
\begin{equation}
z(r,\phi)=\sum_{n=1}^\infty \sqrt{\frac{(4n+2)}{2\pi R^{2n+1}}} r^n\sin n\phi
\end{equation}
However, although this looks all tidy and nice, when I plug $z(r,\phi)$ into $z_{rr}+\frac{1}{r}z_r+\frac{1}{r^2}z_{\phi\phi}=0$, I don't get zero.
Any ideas what has gone wrong?
 A: The Lagrange multiplier condition for extremizing this functional subject to these constraints is $\nabla^2 z = -\lambda z$, as you found. I can derive this from scratch if desired but it seems unnecessary since you found it yourself anyway.
Now the whole problem is radially symmetric, so we assume a radial solution $z(x,y)=Z(r)$. The constraint reads $2 \pi \int_0^R r Z(r)^2 dr = 1$. Our equation reads $Z''+Z'/r=-\lambda Z$, or $r^2 Z'' + r Z' + \lambda r^2 Z = 0$.
First take a detour to note that there is no solution to this problem with $\lambda=0$ (because the equation and BC imply $Z \equiv 0$ which violates the constraint). Then note that if we replace $r$ by $cr$ then the first two terms are invariant while the third term is multiplied by $c^2$. So you can replace $r$ by $s=r/\sqrt{\lambda}$ to get $s^2 Z'' + s Z' + s^2 Z = 0$. Therefore $Z(r)=c_1 J_0(r/\sqrt{\lambda})+c_2 Y_0(r/\sqrt{\lambda})$ where $J_0$ and $Y_0$ are Bessel functions. You can use the boundary conditions to determine $c_1,c_2$ and $\lambda$. (They're not unique.)
If the problem were not radially symmetric then you would need to separate variables.
