Analytic solution to Poisson equation I need to find the analytic solution to this equation, in order to compare it with solution I get from using a numerical solution. However, I have not been able to find the solution. 
I think I can't use separation of variables since the R.H.S cannot be separated. I tried solving it using Green's theorem but to no luck. Any ideas how to approach this?
$\bigtriangleup u = 2x(y −1)(y −2x + x y +2)e^{x−y}$, $(x, y) ∈ (0, 1)×(0, 1)$
Boundary conditions
$u(x, 0) = u(x, 1) = 0$, $x ∈ [0, 1]$
$u(0, y) = u(1, y) = 0$, $y ∈ [0, 1]$
 A: I found the solution:
Since 
$\bigtriangleup u = 2x(y −1)(y −2x + x y +2)e^{x−y}$, $(x, y) ∈ (0, 1)×(0, 1)$
The function u(x,t) must have $e^{x−y}$
From the boundary conditions:


*

*$u(x, 0) = 0$ then u(x,t) must be 0 for x=0, so we assume u(x,t) must
have $e^{x−y}*x$

*$u(x, 1) = 0$ then u(x,t) must be 0 for x=1, so we assume u(x,t) must
have $e^{x−y}*x*(1-x)$

*$u(0, y) = 0$ then u(x,t) must be 0 for y=0, so we assume u(x,t) must
have $e^{x−y}*x*(1-x)*y$

*$u(1, y) = 0$ then u(x,t) must be 0 for y=1, so we assume u(x,t) must
have $e^{x−y}*x*(1-x)*y*(1-y)$


The solution should look like this: $u(x,t)=e^{x−y}*x*(1-x)*y*(1-y)$
Doing the calculations leads to:
$\bigtriangleup u = 2x(y −1)(y −2x + x y +2)e^{x−y}$
Which due to the uniqueness of the solution, is the solution.
Calculations can be seen below:

A: Note that the solution of $-\Delta u = f$, with homogeneous dirichlet BC's is given by $$u(x)=\Phi*f=\int_\Omega\Phi(x-y)f(y)\,dy$$ where $\Phi$ is the fundamental solution of the Laplace equation, and (in 2D) $$\Phi(x)=\frac{\ln(|x|)}{2\pi}.$$
Note that in the sense of distributions, $-\Delta\Phi(x)=\delta(x)$, and so 
$$-\Delta u = \int_\Omega-\Delta\Phi(x-y)f(y)\,dy = \int_\Omega\delta(x-y)f(y)\,dy = f(x).$$
A: with free CAS Maxima:
(%i1)   load(diffops);

(%o1)   "C:\Users\aleksas\maxima\diffops.mac"

(%i2)   eq:'diff(u,x,2)+'diff(u,y,2)=2*x*(y-1)*(y-2*x+x*y+2)*exp(x-y);

(eq)    'diff(u,y,2)+'diff(u,x,2)=2*x*(y-1)*(x*y+y-2*x+2)*%e^(x-y)

(%i3)   partsolpde(eq,u,[x,y]);

(%o3)   x^2*y^2*%e^(x-y)-x*y^2*%e^(x-y)-x^2*y*%e^(x-y)+x*y*%e^(x-y)

(%i4)   factor(%);

(%o4)   (x-1)*x*(y-1)*y*%e^(x-y)

This is solution 
