Solving the Helmholtz equation on a square with Neumann boundary conditions . I would like to solve 
$$\nabla^2u+u=0$$
on the square defined by $[-1,1] \times [-1,1]
 $ subject to $\frac{\partial u}{\partial n} = 1$ on the boundary. By trying solutions of the form $u(x,y) = f(x)g(y)$ I obtain the expression
$$\lambda = \frac{f''-f}{f}= -\frac{g''}{g}$$
Assuming that $\lambda > 0$ I obtain an expression for $g$:
$$g(y) = A\cos\sqrt\lambda(y-1) + B\cos\sqrt\lambda(y+1)$$
where 
$$B=A = -\frac{1}{\sqrt \lambda\sin2\sqrt \lambda}$$
was obtained from the boundary condition $g'(\pm 1) =\pm 1$.
For $f$:
$$f(x) = C\cosh\sqrt{\lambda+1}(y-1) + D\cosh\sqrt{\lambda+1}(y+1)$$
where
$$C=D = \frac{1}{\sqrt {\lambda +1}\sinh2\sqrt {\lambda +1}}$$
by using the boundary conditions  $f(\pm 1) = \pm 1$
I have two questions:


*

*How do I obtain $\lambda$?

*Am I proceeding in the correct way?
 A: I'm assuming $\frac{\partial u}{\partial n}$ refers to the outward-pointing normal (you can reverse the sign if need be), then the boundary conditions, explicitly, are as follows
$$ \begin{array}{ll}
\dfrac{\partial u}{\partial x} = 1, && x = 1 \\
\dfrac{\partial u}{\partial x} = -1, && x = -1 \\
\dfrac{\partial u}{\partial y} = 1, && y = 1 \\
\dfrac{\partial u}{\partial y} = -1, && y = -1 
\end{array} $$
These boundary conditions are non-homogeneous. To solve this, we have to separate the solution into 4 different "pieces"
$$ u(x,y) = u_1(x,y) + u_2(x,y) + u_3(x,y) + u_4(x,y) $$
such that each piece is only non-homogeneous on one boundary and homogeneous on the other 3. Something like this
\begin{array}{ll} \dfrac{\partial u_1}{\partial x} = 1, && x = 1 \\
\dfrac{\partial u_1}{\partial x} = 0, && x = -1 \\
\dfrac{\partial u_1}{\partial y} = 0, && y = \pm 1
\end{array}
and likewise for the other 3 pieces

Here's a basic solution for the first piece. Let $u_1(x,y) = f_1(x)g_1(y)$, the equation separates to
$$ \frac{{f_1}''}{f_1} + 1 = -\frac{{g_1}''}{g_1} = \lambda $$
and the 3 homogeneous B.C.'s force ${f_1}'(-1) = {g_1}'(\pm 1) = 0$
Let's start first with the $y$ part. You can check that only $\lambda \ge 0$ leads to a valid solution. For convenience, we shift the solution to one of the boundaries
$$ g_1(y) = A\cos(\sqrt{\lambda} (y+1)) + B\sin(\sqrt{\lambda}(y+1)) $$
${g_1}'(-1)=0$ obviously gives $B=0$. ${g_1}'(1) = 0$ gives $\sin (2\sqrt{\lambda}) = 0$, therefore $2\sqrt{\lambda} = n\pi$. So we get
$$ {g_1}(y) = \cos \left(\frac{n\pi}{2} (y+1)\right) $$
For the $x$ part, observe that
$$ \lambda - 1 = \frac{n^2\pi^2}{4} - 1 $$
is negative for $n = 0$ and positive otherwise, so we have two cases
$$ f_1(x) = \cos (x+1), \quad n = 0 $$
$$ f_1(x) = \cosh \left(\sqrt{\frac{n^2\pi^2}{4} - 1}(x+1)\right), \quad n = 1, 2, 3, \cdots $$
The law of superposition allows the general solution as
$$ u_1(x,y) = A_0\cos(x+1) + \sum_{n=1}^\infty A_n\cosh\left( \sqrt{\frac{n^2\pi^2}{4}-1}(x+1)\right) \cos\left(\frac{n\pi}{2}(y+1)\right) $$
Applying the non-homogeneous B.C.
$$ \left.\frac{\partial u_1}{\partial x}\right|_{(1,y)} = 1 = -A_0\sin(2) + \sum_{n=1}^\infty A_n \sqrt{\frac{n^2\pi^2}{4}-1}\ \sinh \left( 2\sqrt{\frac{n^2\pi^2}{4}-1}\right)\cos\left(\frac{n\pi}{2}(y+1)\right) $$
This is simply a cosine series in $y$, where the coefficients are
$$ -A_0\sin(2) = \frac{1}{2}\int_{-1}^1 dy = 1 $$
$$ A_n \sqrt{\frac{n^2\pi^2}{4}-1} \sinh \left( 2\sqrt{\frac{n^2\pi^2}{4}-1}\right) = \int_{-1}^1 \cos\left(\frac{n\pi}{2}(y+1)\right) dy = 0 $$
All of the above shows
$$ u_1(x,y) = -\frac{\cos(x+1)}{\sin (2)} $$

Going through the same process for the other 3 pieces, we obtain the final solution
$$ u(x,y) = -\frac{\cos(x+1)+\cos(1-x)+\cos(y+1)+\cos(1-x)}{\sin (2)} $$
A: You are completely off base, but it took me a while to see through this problem as well. The way to start is normally to set the boundary conditions to $0$ at $3$ of the boundaries and then to satisfy the actual boundary condition at the fourth. Here I chose $-\frac{\partial u}{\partial y}=0$ at $y=-1$, $\frac{\partial u}{\partial y}=0$ at $y=1$, $\frac{\partial u}{\partial x}=0$ at $x=1$, but $-\frac{\partial u}{\partial x}=1$ at $x=-1$. Then when you separate variables $u(x,y)=f(x)g(y)$ you get
$$\frac1f\frac{d^2f}{dx^2}=-\lambda=-\frac1g\frac{d^2g}{dy^2}-1$$
And normally you would find the values of $\lambda=1-\frac{n^2\pi^2}4$ that are consistent with the boundary conditions and then find the corresponding $g_n(y)$ and $f_n(x)$ functions and then find the Fourier coefficients by applying the boundary condition at $x=0$. But here already one of the $g_n(y)$ functions is constant so that is the only one that can match the boundary condition. So we have for $1-\lambda=0$
$$g_0^{\prime\prime}(y)=0$$
So
$$g_0(y)=c_1+c_2y$$
$$g_0^{\prime}(-1)=g_0^{\prime}(1)=c_2=0$$
So we have the solution $g_0(y)=1$. Then
$$f_0^{\prime\prime}(x)+f_0(x)=0$$
So
$$f_0(x)=c_1\cos(1-x)+c_2\sin(1-x)$$
$$f_0^{\prime}(1)=-c_2=0$$
So a solution is $u_0(x,y)=a_0f_0(x)g_0(y)=a_0\cos(1-x)$. Then we want
$$-\frac{\partial}{\partial x}u_0(x,y)=-a_0\sin(1-x)=-a_0\sin(2)=1$$
So our solution is
$$u_0(x,y)=-\frac{\cos(1-x)}{\sin(2)}$$
At this point we can see what the solutions will be to satisfy the other $3$ boundary conditions and we have
$$u(x,y)=-\frac{\cos(1-x)}{\sin(2)}-\frac{\cos(1-y)}{\sin(2)}-\frac{\cos(1+x)}{\sin(2)}-\frac{\cos(1+y)}{\sin(2)}$$
So it turned out much easier than anticipated!
