Solving green function for first quadrant with boundary conditions (dirichlet problem?) Let $\Omega = \{x_1,x_2\in\mathbb{R}^2: x_1>0, x_2>0\}$. Solve the problem $$\Delta u = 0 \mbox{ in $\Omega$ }, u\in C^2(\Omega)\cap C(\overline\Omega) \mbox{ bounded}\\u(x_1,0)=u_0(x_1),x_1\ge 0\\u(0,x_2) = 0, x_2\ge 0$$
Where $u_0$ is continuous, bounded from $0$ to $\infty$ and $u_0(0) = 0$
I opened a bounty in Green function of the first quadrant because I thought it would help be solve this problem which is a slighty different one but I don't think it helps. 
I truly have no idea how to modify the answer there to obtain these boundary conditions. 
UPDATE:
The only thing I could think of is to solve a system of dirichlet problems:
$$\Delta u_1 = 0 \mbox{ in } \mathbb{R}_+^n\\ u_1 = -u_2 + u_0 \mbox{ in } \partial \mathbb{R}_+^n$$
$$\Delta u_2 = 0 \mbox{ in } \mathbb{R}_{++}^n\\ u_2 = -u_1\mbox{ in } \partial \mathbb{R}_{++}^n$$
Where $\mathbb{R}_+^n$ is the upper half plane, and $\mathbb{R}_{++}^n$ is the 'upper right plane'. That is $\{(x_1,x_2)| x_1>0\}$
Then if we sum these solutions, we have
$u_3 = u_1 + u_2 \implies \Delta u_3 = 0$ in $\{(x_1,x_2), x_1>0, x_2>0\}$
and at the first border, which is $\{(x_1,0), x_1>0\}$, we have
$u_3 = u_1 + u_2 = -u_2 + u_0 + u_2 = u_0$
and at the second border, which is $\{(0,x_2), x_2>0\}$, we have
$u_3 = u_1 + u_2 = u_1 + -u_1 = 0$
So we should find the solutions of $u_1$ and $u_2$ explictly. I found them to be
$$u_1(x) = \frac{x_2}{\pi}\int \frac{1}{|y-x|^2}(-u_2(x)+u_0(x))\ dy$$
$$u_2(x) = \frac{x_1}{\pi}\int \frac{1}{|y-x|^2}(-u_1(x))\ dy$$
so the solution would be
$$u_1 + u_2 = \\ \frac{x_2}{\pi}\int \frac{1}{|y-x|^2}\left(-\frac{x_1}{\pi}\int \frac{1}{|y-x|^2}(-u_1(x))\ dy + u_0(x)\right)\ dy + \frac{x_1}{\pi}\int \frac{1}{|y-x|^2}(-u_1(x))\ dy$$
Is it right? 

You could use separation of variables and Fourier analysis to obtain solutions. To do this, let $u(x,y)=X(x)Y(y)$, and solve
$$
         -\frac{X''}{X}=\lambda = \frac{Y''}{Y},\\ X(0)=0.
$$
Restricting $Y$ to be bounded gives solutions
$$
            X(x)Y(y)=\sin(\sqrt{\lambda}x)e^{-\sqrt{\lambda}y},\;\;\lambda > 0.
$$
A linear combination of these can be formed using a Fourier integral:
$$
           u(x,y)=\int_{0}^{\infty}c(s)\sin(sx)e^{-sy}ds
$$
In order to specify $u(x,y)=u_0(x)$, the coefficient function $c(s)$ must satisfy
$$
           u_0(x)=\int_{0}^{\infty}c(s)\sin(sx)ds.
$$
Using the Fourier sine transform isolates $c(s)$:
$$
        \frac{2}{\pi}\int_{0}^{\infty}u_0(x)\sin(sx)dx = c(s).
$$
Therefore,
$$
      u(x,y)=\frac{2}{\pi}
\int_{0}^{\infty}\left(\int_{0}^{\infty}\sin(sx')u_0(x')dx'\right)\sin(sx)e^{-sy}ds
$$
