Here's a Poisson equation, $D$ is a unit square i.e. $[0,1]\times[0,1]$. \begin{equation*} \begin{aligned} \Delta u(x,y) &= \varphi(x,y) \quad (x,y)\in D \\ u(x,y) &= 0 \quad (x,y)\in\partial D \\ \end{aligned} \end{equation*}
$\varphi(x,y)=xy(1-x)(1-y)$, the question asked me to solve it with separate variable method. I truly did know how to solve it when $\Delta u=0$, but I don't know how to use separate variable method to solve this inhomogeneous one.
I searched for some similar questions, maybe it's a better way to use charasteristic function,
$\varphi(x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{mn}\sin{m\pi x}\sin{n\pi y}$
I would like to know how did this method worked in the problem and why it can be expressed like that.
Firstly consider the homogeneous Laplacian equation $\Delta u=0$. Use separate variable method, $u(x,y)=X(x)Y(y)$, then $\frac{X''}{X}=-\frac{Y''}{Y}=c$. Obviously $c$ is negative, $X(x)=A\sin\sqrt{c}x+B\cos\sqrt{c}x$ where $A,B$ are constants. $X(0)=0\Rightarrow B=0$, $X(1)=0\Rightarrow c=n^2\pi^2(n\in\mathbb{N})$, $X(x)=A\sin{n\pi x}$. Turn to $Y(y)$ to get the similar result (but needn't to plug $c=n^2\pi^2$ in), $Y(y)=C\sin{m\pi y} (m\in\mathbb{N})$. This is why $u(x,y)$ can have the form $\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}A_{nm}\sin{n\pi x}\sin{m\pi y}$.
In the inhomogeneous situation, it is time to use separate method to $\varphi(x,y)$. In PDE's theory a solution to Poisson equation can have the form $\varphi(x,y)=\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}b_{nm}\sin{n\pi x}\sin{m\pi y}$, here $b_{nm}$ is the Fourier coefficient similar to the homogeneous one.
It's easy to know when make a Fourier expansion to $\varphi(x,y)$ with respect to $x$, it is $\varphi(x,y)=\sum_{n=1}^{\infty}a_{n}(y)\sin{n\pi x}$, here $a_{n}(y)$ is \begin{equation} a_{n}(y)=2\int_{0}^{1}\varphi(x,y)\sin{n\pi x}\ \mathrm{d}x \end{equation}
Similarly do the same with respect to $y$ above, $a_{n}(y)=\sum_{m=1}^{\infty}b_{nm}\sin{m\pi y}$.
Calculate that to get \begin{equation} b_{nm}=4\int_{0}^{1}\left(\int_{0}^{1}\varphi(x,y)\sin{m\pi y}\ \mathrm{d}y\right)\sin{n\pi x}\ \mathrm{d}x \end{equation}
Note that $\Delta u=-(m^2+n^2)\pi^2 u=\varphi$, so it's easy to get $A_{nm}$: \begin{equation} A_{nm}=\frac{-b_{nm}}{(m^2+n^2)\pi^2} \end{equation}
Finally we have the solution \begin{equation} u(x,y)=\frac{-4}{\pi^2}\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\frac{1}{m^2+n^2}\left(\int_{0}^{1}\left(\int_{0}^{1}\varphi(x,y)\sin{m\pi y}\ \mathrm{d}y\right)\sin{n\pi x}\ \mathrm{d}x\right)\sin{n\pi x}\sin{m\pi y} \end{equation}
Note that $\varphi(x,y)=xy(1-x)(1-y)$ is already separated, the solution can be more simple: \begin{equation} u(x,y)=\frac{-4}{\pi^2}\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\frac{1}{m^2+n^2}\left(\int_{0}^{1}x(1-x)\sin{n\pi x}\ \mathrm{d}x\right)\left(\int_{0}^{1}y(1-y)\sin{m\pi y}\ \mathrm{d}y\right) \end{equation}