How do I solve an inhomogeneous Helmholtz boundary value problem in 2-D with a rectangular boundary?

I need to solve the following BVP:

$\Delta u - 1/\delta * u = R(x,y)$,

where $\Delta$ is the laplacian operator. and boundary conditions:

$u(0,y)=u(L,y)=u(x,0)=u(x,L)=0$

where $L=1$, $\delta=0.01$, and $R(x,y)=1$.

I understand the basic concept of separation of variables, i.e. rewriting $u(x,y)$ as the product of some function of $x, X(x)$, and some function of $y, Y(y)$. Then by plugging $u(x,y)=X(x)Y(y)$ into the PDE, obtaining two ODEs that are much easier to solve. My solution started as...

$X''Y+XY''-1/\delta*XY=0$

$X''/X=(1/\delta*Y-Y'')/Y=-\lambda^2$

$X(x) = Asin\lambda x + Bcos\lambda y$

$Y(y) = C_1 e^{\sqrt{1/\delta+\lambda^2}y} + C_2 e^{-\sqrt{1/\delta+\lambda^2}y}$

eventually...

$u(x,y) = \sum{A_n(e^{\sqrt{1/\delta+\lambda^2}y} + e^{-\sqrt{1/\delta+\lambda^2}y})sin{(n \pi x)}}$

However, this was merely for the homogeneous case, not the BVP as given. And in addition, it seemed counter intuitive to me that the functions $X(x)$ and $Y(y)$ should be of different forms. Why would $X(x)$ be composed of sinusoidal functions and $Y(y)$ be compased of exponential terms if $u(x,y)$ is taken to the second partial derivative in both the x and y directions?

I am a little suspicious of your solution to the homogeneous problem, since $u(x, y) \equiv 0$ is a solution to this problem too (with the boundary conditions) --- is there a reason in particular that this is not the unique solution to your problem? Have you actually computed the $A_n$?

As for your harder question: how to solve the inhomogeneous problem: this requires the use of Green's functions as well as the use of a Sommerfeld radiation condition As seen here. This is almost too much to explain here, since Green's functions usually occupy at least a quarter, if not three of applied math graduate school.

A good book for practical use of Green's function is Stakgold - this is what I learned from.

My application of separation of variables was incorrect:

$\Delta u - 1/\delta * u = R(x,y)$,

Consider the homogeneous case:

$X''Y+XY''-1/\delta*XY=0$

$\frac{X''Y}{XY}+\frac{XY''}{XY}-1/\delta*\frac{XY}{XY}=0$

$\frac{X''}{X}+\frac{Y''}{Y}+1/\delta=0$

Because X is only a function of x and Y is only a function of y, the only way for the above to be true is if $\frac{X''}{X}$ and $\frac{Y''}{Y}$ are constants.

$\frac{X''}{X}=-\lambda^2$

$\frac{Y''}{Y}=-\gamma^2$

$X(x) = Asin\lambda x + Bcos\lambda y$

$X(x) = Csin\gamma x + Dcos\gamma y$

eventually, given the boundary conditions, the general solution is...

$u(x,y) = \sum{A_n sin(n\pi x) sin(n\pi y)}$