How to solve the following Poisson equation? Solve the following Poisson equation $$u_{xx}+u_{yy}=x^2+4xy  \text{ for }  x^2+y^2<4\\
u(x,y)=-2 \text{ for } x^2+y^2=4$$
I know how to solve homogenous Laplace equation but I have no idea how to solve can you help
 A: You can just set up a discretized linear equation system.
Make ansatz of truncated power series expansion $$u(x,y) = \sum_{\forall i,j \in\{0,\cdots,N\}} c_{ij}x^iy^j$$
Vectorize $c_{ij}$ lexicographically. Implement $\bf D$ as $\{0,1,\cdots N-1\}$ 1 step off-diagonal matrix. Example for 6th degree polynomial :
$${\bf D} = \left[\begin{array}{lllllll}&1&&&&&\\&&2&&&&\\&&&3&&&\\&&&&4&&\\&&&&&5&\\&&&&&&6\\&&&&&&\end{array}\right]$$
Now we can construct 2D partial differentiators like so:
$${\bf D_x = ({\bf I}\otimes {\bf D})}\\{\bf D_y = (D\otimes I)}$$
Now set up boundary conditions. We can sample uniformly $$x+yi = 2\exp(2\pi i k/N), \forall k \in\{0,\cdots,N-1\}$$
Remains to do is to build equation system and solve it. Here is what a solution can look like (the blue circles are the linearly spaced points used to ensure boundary conditions are followed) :

Just for verification purposes now that we have an exact solution by @Yalikesifulei here is a truncated representation of the numerical solution (any number abs val $<10^{-10}$ has been set to $0$) :
$$\left[\begin{array}{lllll}-2.500&0&0.166667&0&-0.01041666\\0&-1.3333&0&0.333333&0\\-0.166667&0&0.06249999&0&0\\0&0.33333&0&0&0\\0.07291666&0&0&0&0\end{array}\right]$$
The three trickiest coefficients which we may not see straight from ocular inspection are the same are
for $y^4$ :
$1/12-1/96 = 0.07916666...$
and $x^4$ :
$-1/96 = -0.010416666...$
and for $x^2y^2$ :
$1/12-2/96 = 0.0625$
A: It is beneficial to adopt the plane polar coordinates for this problem. One has
\begin{align}
\frac{1}{r}\frac{\partial}{\partial{r}}\Bigg(r\frac{\partial}{\partial{r}}\psi\Bigg) +\frac{1}{r^{2}}\frac{\partial^{2}}{\partial{\varphi}^{2}}\psi=r^{2}[\cos^{2}\varphi+2\sin2\varphi], \quad \psi(2, \varphi)=-2.
\end{align}
The solution to the homogenous equation obeys
\begin{align}
\frac{1}{r}\frac{\partial}{\partial{r}}\Bigg(r\frac{\partial}{\partial{r}}\psi_{h}\Bigg) +\frac{1}{r^{2}}\frac{\partial^{2}}{\partial{\varphi}^{2}}\psi_{h}=0
\end{align}
We let
\begin{align}
\psi_{h}(r, \varphi)=\frac{1}{\sqrt{2\pi}}\sum_{n\in\mathbb{Z}}\hat\psi_{h}^{(n)}(r)e^{in\varphi},
\end{align}
so that
\begin{align}
r^{2}\frac{\partial^{2}}{\partial{r}^{2}}\hat\psi_{h}^{(n)}+r\frac{\partial}{\partial{r}}\hat\psi_{h}^{(n)}-n^{2}\hat\psi_{h}^{(n)}=0.
\end{align}
This is a standard Euler equation which has the solution
\begin{align}
\hat\psi_{h}^{(n)}(r)=C_{h, +}^{(n)}r^{n}+C_{h, -}^{(n)}r^{-n}.
\end{align}
Hence
\begin{align}
\psi_{h}(r, \varphi)=\sum_{n\in\mathbb{Z}}[C_{h, +}^{(n)}r^{n}+C_{h, -}^{(n)}r^{-n}]\frac{e^{in\varphi}}{\sqrt{2\pi}}.
\end{align}
Now we consider the inhomogeneous part, i.e. the particular solution
\begin{align}
\frac{1}{r}\frac{\partial}{\partial{r}}\Bigg(r\frac{\partial}{\partial{r}}\psi_{p}\Bigg) +\frac{1}{r^{2}}\frac{\partial^{2}}{\partial{\varphi}^{2}}\psi_{p}=r^{2}[\cos^{2}\varphi+2\sin2\varphi], \quad \psi(2, \varphi)=-2.
\end{align}
Again we Fourier expand
\begin{align}
\frac{1}{r}\frac{\partial}{\partial{r}}\Bigg(r\frac{\partial}{\partial{r}}\psi_{p}^{(n)}\Bigg) -\frac{n^{2}}{r^{2}}\psi_{p}^{(n)}=r^{2}\int_{-\pi}^{\pi}d\varphi[\cos^{2}\varphi+2\sin2\varphi]\frac{e^{-in\varphi}}{\sqrt{2\pi}}.
\end{align}
We have
\begin{align}
\int_{-\pi}^{\pi}d\varphi[\cos^{2}\varphi+2\sin2\varphi]\frac{e^{-in\varphi}}{\sqrt{2\pi}}=\frac{\sqrt{2\pi}}{4}[2\delta_{n, 0}+[1-4i]\delta_{n, 2}+[1+4i]\delta_{n, -2}].
\end{align}
Hence the "non-trivial" modes satisfy
\begin{align}
\frac{\partial^{2}}{\partial{r}^{2}}\psi_{p}^{(0)}+\frac{1}{r}\frac{\partial}{\partial{r}}\psi_{p}^{(0)}=&r^{2}\frac{\sqrt{2\pi}}{2},\\
\frac{\partial^{2}}{\partial{r}^{2}}\psi_{p}^{(2)}+\frac{1}{r}\frac{\partial}{\partial{r}}\psi_{p}^{(2)} -\frac{4}{r^{2}}\psi_{p}^{(2)}=&r^{2}\frac{[1-4i]\sqrt{2\pi}}{4},\\
\frac{\partial^{2}}{\partial{r}^{2}}\psi_{p}^{(-2)}+\frac{1}{r}\frac{\partial}{\partial{r}}\psi_{p}^{(-2)} -\frac{4}{r^{2}}\psi_{p}^{(-2)}=&r^{2}\frac{[1+4i]\sqrt{2\pi}}{4}.
\end{align}
We let
\begin{align}
\psi_{p}^{(0)}=&\frac{\sqrt{2\pi}}{32}r^{4},\\
\psi_{p}^{(2)}=&\frac{[1-4i]\sqrt{2\pi}}{48}r^{4},\\
\psi_{p}^{(-2)}=&\frac{[1+4i]\sqrt{2\pi}}{48}r^{4}.
\end{align}
It follows
\begin{align}
\psi_{p}(r, \varphi)=&\frac{r^{4}}{8}\Bigg[\frac{1}{4}+\frac{2}{3}\cos(2\varphi)+\frac{8}{3}\sin(2\varphi)\Bigg].
\end{align}
The general solution is thus given by
\begin{align}
\psi(r, \varphi)=\sum_{n\in\mathbb{Z}}[C_{h, +}^{(n)}r^{n}+C_{h, -}^{(n)}r^{-n}]\frac{e^{in\varphi}}{\sqrt{2\pi}}+\psi_{p}(r, \varphi).
\end{align}
The boundary conditions dictate
\begin{align}
\sum_{n\in\mathbb{Z}}[C_{h, +}^{(n)}2^{n}+C_{h, -}^{(n)}2^{-n}]\frac{e^{in\varphi}}{\sqrt{2\pi}}+\psi_{p}(2, \varphi)=-2.
\end{align}
We multiply by $\frac{e^{-in\varphi}}{\sqrt{2\pi}}$ an integrate over $2\pi$, so that
\begin{align}
C_{h, +}^{(n)}2^{n}+C_{h, -}^{(n)}2^{-n}+\psi_{p}^{(0)}(2)\delta_{n, 0}+\psi_{p}^{(2)}(2, \varphi)\delta_{n, 2}+\psi_{p}^{(-2)}(2)\delta_{n, -2}=-2\sqrt{2\pi}\delta_{n, 0}.
\end{align}
We have in particular
\begin{align}
C_{h, +}^{(0)}+C_{h, -}^{(0)}=&-\frac{5}{2}\sqrt{2\pi},\\
C_{h, +}^{(2)}4+C_{h, -}^{(2)}\frac{1}{4}=&-\frac{[1-4i]\sqrt{2\pi}}{3},\\
C_{h, +}^{(-2)}\frac{1}{4}+C_{h, -}^{(2)}4=&-\frac{[1+4i]\sqrt{2\pi}}{3},\\
C_{h, +}^{(n)}=-C_{h, -}^{(n)}2^{-2n},\quad n\neq0, \ \pm2.
\end{align}
This completes the solution.
A: The other answers are good. I'd like to give my own thoughts.
You can approach this problem in a similar way as I do here, by breaking the problem into simpler parts. Your problem is
$$\begin{cases}
(\Delta u)( x,y) =x^{2} +4xy & ( x,y) \in \mathbb{B}( 0,4)\\
u( x,y) =-2 & ( x,y) \in \partial \mathbb{B}( 0,4)
\end{cases}\tag{0}$$
Where $\Delta=\partial_x^2+\partial_y^2$.
I would suggest first finding a solution (call it $u_1$) to
$$\begin{cases}
(\Delta u)( x,y) =x^{2} +4xy & ( x,y) \in \mathbb{B}( 0,4)\\
u( x,y) =0 & ( x,y) \in \partial \mathbb{B}( 0,4)
\end{cases}\tag{1}$$
And then find a solution (call it $u_2$) to
$$\begin{cases}
(\Delta u)( x,y) =0 & ( x,y) \in \mathbb{B}( 0,4)\\
u( x,y) =-2 & ( x,y) \in \partial \mathbb{B}( 0,4)
\end{cases}\tag{2}$$
And then by the linearity of the Laplacian, $u_1+u_2$ will be a solution to $\boldsymbol{(0)}$. What is nice about this approach is that the solution to $\boldsymbol{(2)}$ is already well known. What we can do is make a change of coordinates
$$x'=x/4~~;~~y'=y/4$$
To move the problem to the unit disk, and then use the Poisson Kernel. Switching to polar coordinates, given the BVP on the unit disk
$$\begin{cases}
(\Delta u)(r,\theta) =0 & r<1\\
u(1,\theta)=h(\theta)& \\
(r,\theta)\in[0,1]\times(-\pi,\pi]
\end{cases}\tag{0}$$
The solution is given by
$$u(r,\theta)=\frac{1}{2\pi}\int_{-\pi}^{\pi}h(\phi)\frac{1-r^2}{1+r^2+2r\cos(\theta-\phi)}\mathrm d\phi$$
See pages 160-166 in Peter J Olver's Introduction to Partial Differential Equations for a derivation of this formula.
So in your case you can pick $h=-2$. Unfortunately I don't think $\boldsymbol{(1)}$ is quite this easy, since $x^2+4xy$ doesn't transform all that nicely under polar coordinates. So perhaps a closed form for this part is not possible. I'm sure the other answers will be plenty useful for this bit.
