Added March 13, 2019.
Solve for \(u(x,y) \) \begin {align*} \frac {u_{xx}}{A} + \frac {u_{xx}}{B} = -2 \theta \end {align*}
Where \(A,B,\theta \) are constants, and the boundary conditions are \begin {align*} u(x, -b) &= 0\\ u(x, b) &= 0\\ u(-a, y) &= 0\\ u(a, y) &= 0 \end {align*}
Mathematica ✗
ClearAll["Global`*"]; pde = D[u[x, y], {x, 2}]/A + D[u[x, y], {y, 2}]/B == -2*theta; bc = {u[x, -b] == 0, u[x, b] == 0, u[-a, y] == 0, u[a, y] == 0}; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bc}, u[x, y], {x, y}], 60*10]];
Failed
Maple ✗
restart; pde := diff(u(x,y),x$2)/A+diff(u(x,y),y$2)/B = -2*theta; bc := u(x,-b)=0, u(x,b)=0, u(-a,y)=0, u(a,y) = 0; cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,bc],u(x,y))),output='realtime'));
sol=()
Hand solution
solve \begin {align*} \frac {u_{xx}}{A}+\frac {u_{yy}}{B} & =-2\theta \\ Bu_{xx}+Au_{yy} & =-2\theta AB\\ & =C \end {align*}
Where \(C=-2\theta AB\) is a new constant. With boundary conditions\begin {align*} u\left ( x,-b\right ) & =0\\ u\left ( x,b\right ) & =0\\ u\left ( -a,y\right ) & =0\\ u\left ( a,y\right ) & =0 \end {align*}
To simplify solution, shift the rectangle so its lower left corner on the origin. Let \(\tilde {x}=x+a\), and \(\tilde {y}=y+b\). The boundary conditions becomes\begin {align*} u\left ( \tilde {x},0\right ) & =0\\ u\left ( \tilde {x},2b\right ) & =0\\ u\left ( 0,\tilde {y}\right ) & =0\\ u\left ( 2a,\tilde {y}\right ) & =0 \end {align*}
And the pde becomes \(Bu_{\tilde {x}\tilde {x}}+Au_{\tilde {y}\tilde {y}}=C\). Instead of keep writing \(\tilde {x},\tilde {y}\), will use \(x,y\), but remember that these are shifted version. At the end, we shift back.
Hence the PDE to solve is \(Bu_{xx}+Au_{yy}=C\) with BC\begin {align*} u\left ( x,0\right ) & =0\\ u\left ( x,2b\right ) & =0\\ u\left ( 0,y\right ) & =0\\ u\left ( 2a,y\right ) & =0 \end {align*}
Using eigenfunction expansion method. Let \begin {equation} u\left ( x,y\right ) =\sum _{n=1}^{\infty }b_{n}\left ( y\right ) X_{n}\left ( x\right ) \tag {1} \end {equation} Where \(X_{n}\left ( x\right ) \) is eigenfunctions for \(X^{\prime \prime }\left ( x\right ) +\lambda _{n}X\left ( x\right ) =0\) with boundary conditions \(X\left ( 0\right ) =X\left ( 2a\right ) =0\). This has eigenfunctions as \(X_{n}\left ( x\right ) =\sin \left ( \sqrt {\lambda _{n}}x\right ) \) with eigenvalues \(\lambda _{n}=\left ( \frac {n\pi }{2a}\right ) ^{2}\) for \(n=1,2,\cdots \).
Substituting (1) into the PDE \(Bu_{xx}+Au_{yy}=C\) gives\[ B\sum _{n=1}^{\infty }b_{n}\left ( y\right ) X_{n}^{\prime \prime }\left ( x\right ) +A\sum _{n=1}^{\infty }b_{n}^{\prime \prime }\left ( y\right ) X_{n}\left ( x\right ) =C \] Expanding \(C\) (a constant) as Fourier sine series the above becomes\[ B\sum _{n=1}^{\infty }b_{n}\left ( y\right ) X_{n}^{\prime \prime }\left ( x\right ) +A\sum _{n=1}^{\infty }b_{n}^{\prime \prime }\left ( y\right ) X_{n}\left ( x\right ) =\sum _{n=1}^{\infty }q_{n}X_{n}\left ( x\right ) \] But \(X_{n}^{\prime \prime }\left ( x\right ) =-\lambda _{n}X_{n}\left ( x\right ) \), hence the above becomes\begin {align} -B\sum _{n=1}^{\infty }\lambda _{n}b_{n}\left ( y\right ) X_{n}\left ( x\right ) +A\sum _{n=1}^{\infty }b_{n}^{\prime \prime }\left ( y\right ) X_{n}\left ( x\right ) & =\sum _{n=1}^{\infty }q_{n}X_{n}\left ( x\right ) \nonumber \\ Ab_{n}^{\prime \prime }\left ( y\right ) -B\lambda _{n}b_{n}\left ( y\right ) & =q_{n}\tag {1A} \end {align}
But \begin {align*} C & =\sum _{n=1}^{\infty }q_{n}X_{n}\left ( x\right ) \\ \int _{0}^{2a}CX_{n}\left ( x\right ) dx & =q_{n}\int _{0}^{2a}X_{n}^{2}\left ( x\right ) dx\\ \int _{0}^{2a}C\sin \left ( \sqrt {\lambda _{n}}x\right ) dx & =q_{n}\int _{0}^{2a}\sin ^{2}\left ( \sqrt {\lambda _{n}}x\right ) dx\\ \frac {-C}{\sqrt {\lambda _{n}}}\left ( \left ( -1\right ) ^{n}-1\right ) & =q_{n}a\\ q_{n} & =\frac {-C}{a\sqrt {\lambda _{n}}}\left ( \left ( -1\right ) ^{n}-1\right ) \end {align*}
Hence (1A) becomes\[ Ab_{n}^{\prime \prime }\left ( y\right ) -B\lambda _{n}b_{n}\left ( y\right ) =\frac {-C}{a\sqrt {\lambda _{n}}}\left ( \left ( -1\right ) ^{n}-1\right ) \] This is standard second order linear ODE. The solution is\[ b_{n}\left ( y\right ) =D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}y}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}y}+\frac {C}{aB\lambda _{n}^{\frac {3}{2}}}\left ( \left ( -1\right ) ^{n}-1\right ) \] Using the above in (1) gives the solution\begin {equation} u\left ( x,y\right ) =\sum _{n=1}^{\infty }\left ( D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}y}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}y}+\frac {C}{aB\lambda _{n}^{\frac {3}{2}}}\left ( \left ( -1\right ) ^{n}-1\right ) \right ) X_{n}\left ( x\right ) \tag {1A} \end {equation} We now need to find \(D_{n},E_{n}\).
Case \(n\) even
When \(n\) is even \(\left ( \left ( -1\right ) ^{n}-1\right ) =0\) and the solution (1A) becomes\[ u\left ( x,y\right ) =\sum _{n=1}^{\infty }\left ( D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}y}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}y}\right ) X_{n}\left ( x\right ) \] At \(y=0\) the above gives\[ 0=\sum _{n=1}^{\infty }\left ( D_{n}+E_{n}\right ) \sin \left ( \sqrt {\lambda _{n}}x\right ) \] Therefore\begin {equation} D_{n}+E_{n}=0\tag {2} \end {equation} And at \(y=2b\)\[ 0=\sum _{n=1}^{\infty }\left ( D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}2b}\right ) \sin \left ( \sqrt {\lambda _{n}}x\right ) \] Therefore\begin {equation} D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}2b}=0\tag {3} \end {equation} From (2,3) we see that \(D_{n}=E_{n}=0\), Hence \(u\left ( x,y\right ) =0\) when \(n\) even.
Case \(n\) odd
When \(n\) is odd \(\left ( \left ( -1\right ) ^{n}-1\right ) =-2\) and the solution (1A) becomes\[ u\left ( x,y\right ) =\sum _{n=1}^{\infty }\left ( D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}y}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}y}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\right ) X_{n}\left ( x\right ) \] At \(y=0\) the above gives\[ 0=\sum _{n=1}^{\infty }\left ( D_{n}+E_{n}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\right ) \sin \left ( \sqrt {\lambda _{n}}x\right ) \] Therefore\begin {equation} D_{n}+E_{n}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}=0\tag {4} \end {equation} And at \(y=2b\)\[ 0=\sum _{n=1}^{\infty }\left ( D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}2b}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\right ) \sin \left ( \sqrt {\lambda _{n}}x\right ) \] Therefore\begin {equation} D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}2b}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}=0\tag {5} \end {equation} Solving (4,5) for \(D_{n},E_{n}\) gives\begin {align*} D_{n} & =\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\frac {1}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\\ E_{n} & =\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\frac {e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}} \end {align*}
Therefore the final solution from (1A) becomes\begin {align*} u\left ( x,y\right ) & =\sum _{n=1,3,5,\cdots }^{\infty }\left ( D_{n}e^{\sqrt {\frac {B}{A}\lambda _{n}}y}+E_{n}e^{-\sqrt {\frac {B}{A}\lambda _{n}}y}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\right ) X_{n}\left ( x\right ) \\ & =\sum _{n=1,3,5,\cdots }^{\infty }\left ( \left ( \frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\frac {1}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{\sqrt {\frac {B}{A}\lambda _{n}}y}+\left ( \frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\frac {e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{-\sqrt {\frac {B}{A}\lambda _{n}}y}-\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\right ) \sin \left ( \sqrt {\lambda _{n}}x\right ) \end {align*}
Where \(\lambda _{n}=\left ( \frac {n\pi }{2a}\right ) ^{2}\). Switching back to original coordinates using \(\tilde {x}=x+a\), and \(\tilde {y}=y+b\), then the above is\[ u\left ( x,y\right ) =\sum _{n=1,3,5,\cdots }^{\infty }\left ( \left ( \frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\frac {1}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{\sqrt {\frac {B}{A}\lambda _{n}}\left ( y+b\right ) }+\left ( \frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\frac {e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e\left ( ^{-\sqrt {\frac {B}{A}\lambda _{n}}y+b}\right ) -\frac {2C}{aB\lambda _{n}^{\frac {3}{2}}}\right ) \sin \left ( \sqrt {\lambda _{n}}\left ( x+a\right ) \right ) \] Where \(C=-2\theta AB\), hence\begin {align*} u\left ( x,y\right ) & =\sum _{n=1,3,5,\cdots }^{\infty }\left ( \left ( \frac {-4\theta AB}{aB\lambda _{n}^{\frac {3}{2}}}\frac {1}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{\sqrt {\frac {B}{A}\lambda _{n}}\left ( y+b\right ) }+\left ( \frac {-4\theta AB}{aB\lambda _{n}^{\frac {3}{2}}}\frac {e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{-\sqrt {\frac {B}{A}\lambda _{n}}\left ( y+b\right ) }+\frac {4\theta AB}{aB\lambda _{n}^{\frac {3}{2}}}\right ) \sin \left ( \sqrt {\lambda _{n}}\left ( x+a\right ) \right ) \\ & =\sum _{n=1,3,5,\cdots }^{\infty }\left ( \left ( \frac {-4\theta A}{a\lambda _{n}^{\frac {3}{2}}}\frac {1}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{\sqrt {\frac {B}{A}\lambda _{n}}\left ( y+b\right ) }+\left ( \frac {-4\theta A}{a\lambda _{n}^{\frac {3}{2}}}\frac {e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}{1+e^{\sqrt {\frac {B}{A}\lambda _{n}}2b}}\right ) e^{-\sqrt {\frac {B}{A}\lambda _{n}}\left ( y+b\right ) }+\frac {4\theta A}{a\lambda _{n}^{\frac {3}{2}}}\right ) \sin \left ( \sqrt {\lambda _{n}}\left ( x+a\right ) \right ) \end {align*}
____________________________________________________________________________________
Taken from Mathematica DSolve help pages.
Solve for \(u\left ( x,y\right ) \) \begin {align*} u_{xx}+ u_{yy} & = 6x - 6y \end {align*}
Boundary conditions \begin {align*} u(x, 0) &= 1 + 11 x + x^3\\ u(x, 2) &= -7 + 11 x + x^3\\ u(0, y) &= 1 - y^3\\ u(4, y) &= 109 - y^3 \end {align*}
Mathematica ✓
ClearAll["Global`*"]; pde = Laplacian[u[x, y], {x, y}] == 6*x - 6*y; bc = {u[x, 0] == 1 + 11*x + x^3, u[x, 2] == -7 + 11*x + x^3, u[0, y] == 1 - y^3, u[4, y] == 109 - y^3}; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bc}, u[x, y], {x, y}], 60*10]];
\[\left \{\left \{u(x,y)\to x^3+11 x-y^3+1\right \}\right \}\]
Maple ✓
restart; pde := diff(u(x,y),x$2)+diff(u(x,y),y$2)=6*x-6*y; bc := u(x,0)=1+11*x+x^3, u(x,2)=-7+11*x+x^3, u(0,y)=1-y^3, u(4,y)=109-y^3; cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,bc],u(x,y))),output='realtime'));
\[u \left ( x,y \right ) ={x}^{3}-{y}^{3}+11\,x+1\]
____________________________________________________________________________________
Added January 13, 2020
Solve Poisson PDE \[ \frac {\partial ^{2}f}{\partial x^{2}}+\frac {\partial ^{2}f}{\partial y^{2}}=6y \]
Mathematica ✓
ClearAll["Global`*"]; pde = Laplacian[u[x, y], {x, y}] == 6*y; sol = AbsoluteTiming[TimeConstrained[DSolve[pde, u[x, y], {x, y}], 60*10]];
\[\left \{\left \{u(x,y)\to c_1(x-i y)+c_2(x+i y)+3 x^2 y\right \}\right \}\]
Maple ✓
restart; pde:=VectorCalculus:-Laplacian(u(x,y),'cartesian'[x,y])=6*y; cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve(pde,u(x,y))),output='realtime'));
\[u \left ( x,y \right ) ={\it \_F1} \left ( y-ix \right ) +{\it \_F2} \left ( y+ix \right ) +3\,{x}^{2}y\]
Hand solution
Solve Poisson PDE \[ \frac {\partial ^{2}f}{\partial x^{2}}+\frac {\partial ^{2}f}{\partial y^{2}}=6y \]
The solution is \begin {equation} f=f_{h}+f_{p}\tag {1} \end {equation} Where \(f_{h}\) is the homogenous solution to Laplace \(\frac {\partial ^{2}f}{\partial x^{2}}+\frac {\partial ^{2}f}{\partial y^{2}}=0\) and \(f_{p}\) is a particular solution. The homogeneous solution is easily found as follows. Let \[ f_{h}=F\left ( mx+y\right ) \] Then \(\frac {\partial f_{h}}{\partial x}=mF^{\prime }\) and \(\frac {\partial ^{2}f_{h}}{\partial x^{2}}=m^{2}F^{\prime \prime }\) and \(\frac {\partial f_{h}}{\partial y}=F^{\prime }\) and \(\frac {\partial ^{2}f_{h}}{\partial y^{2}}=F^{\prime \prime }\). Substituting these into \(\frac {\partial ^{2}f_{h}}{\partial x^{2}}+\frac {\partial ^{2}f_{h}}{\partial y^{2}}=0\) gives\begin {align*} m^{2}F^{\prime \prime }+F^{\prime \prime } & =0\\ m^{2}+1 & =0\\ m & =\pm i \end {align*}
Since we assumed \(f_{h}=F\left ( mx+y\right ) \), then the homogeneous solution is sum of two arbitrary functions (one for each root of \(m\))\[ f_{h}=F_{1}\left ( ix+y\right ) +F_{2}\left ( -ix+y\right ) \] To find particular solution, let \begin {equation} f_{p}=Ax^{n}y^{m}\tag {2} \end {equation} Hence \(\frac {\partial f_{p}}{\partial x}=Anx^{n-1}y^{m}\) and \(\frac {\partial ^{2}f_{p}}{\partial x^{2}}=An\left ( n-1\right ) x^{n-2}y^{m}\) and \(\frac {\partial f_{p}}{\partial y}=Ax^{n}my^{m-1}\) and \(\frac {\partial ^{2}f_{p}}{\partial y^{2}}=Ax^{n}m\left ( m-1\right ) y^{m-2}\). Substituting these into \(\frac {\partial ^{2}f_{p}}{\partial x^{2}}+\frac {\partial ^{2}f_{p}}{\partial y^{2}}=6y\) gives\begin {align*} An\left ( n-1\right ) x^{n-2}y^{m}+Ax^{n}m\left ( m-1\right ) y^{m-2} & =6y\\ y^{m}\left ( An\left ( n-1\right ) x^{n-2}+Ax^{n}m\left ( m-1\right ) y^{-2}\right ) & =6y \end {align*}
Comparing terms, then \(m=1\) and \[ An\left ( n-1\right ) x^{n-2}+Ax^{n}m\left ( m-1\right ) y^{-2}=6 \] Since \(m=1\) then the above simplifies to\[ An\left ( n-1\right ) x^{n-2}=6 \] Since there is no \(x\) in RHS, then \(n-2=0\) or \(n=2\) and the above becomes\begin {align*} 2A\left ( 2-1\right ) & =6\\ A & =3 \end {align*}
Hence (2) becomes\[ f_{p}=Ax^{n}y^{m}=3x^{2}y \] And the complete solution (1) is\begin {align*} f\left ( x,y\right ) & =f_{h}+f_{p}\\ & =F_{1}\left ( ix+y\right ) +F_{2}\left ( -ix+y\right ) +3x^{2}y \end {align*}