(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}
\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?
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.