A specific case of Laplace's equation I'm quite new in the field of differential equations.
Could anyone recommend, how can I solve Laplace's equation of the following form:
\begin{alignat}{4}
\Delta u (x,y) &= f(u(x_0, y_0)) &\quad (x,y) &\in &\Omega, \\
u(x,y) &= g(x,y) &\quad (x,y) &\in \partial &\Omega,
\end{alignat}
where $f$ is a nonlinear function and $(x_0, y_0) \in \Omega$ is a fixed point. 
What bothers me is that right-hand side depends on the value of the unknown function $u$ at the specific point $(x_0,y_0)$.
How can one solve such equations analytically or numerically (e.g. using the finite-element method)?
 A: In reality, your problem consists of solving a single nonlinear equation of the form $$K(c) = 0$$ with respect to $c \in \mathbb{R}$. The exact form of $K$ is given below. In principle, the solution can be accomplished using the robust secant method, but locating the initial bracket might be difficult.
Now, let $f : I \rightarrow \mathbb{R}$ and let $w: I \times \Omega \rightarrow \mathbb{R}$ be given by $$w(c,x,y) = u(x,y),$$ where $u : \Omega \rightarrow \mathbb{R}$ is the unique solution of Laplace's equation 
\begin{alignat}{4} 
\Delta u(x,y) &= f(c), &\quad (x,y) &\in &\Omega, \\
u(x,y) &= g(x,y), &\quad (x,y) &\in \partial &\Omega.
\end{alignat}
The critical function $K : I \rightarrow \mathbb{R}$ is given by $$K(c) = w(c,x_0,y_0) - c.$$ This equation may or may not have a solution, but assume that $K(\xi) = 0$ and let $u$ denote the solution of
\begin{alignat}{4} 
\Delta u(x,y) &= f(\xi), &\quad (x,y) &\in &\Omega, \\
u(x,y) &= g(x,y), &\quad (x,y) &\in \partial &\Omega.
\end{alignat}
Then you are done, because $$f(\xi) = f(w(\xi,x_0,y_0)) = f(u(x_0,y_0)).$$
In practice, there is the added problem of solving the equation $K(c) = 0$ in a reasonable amount of time. To that end, I propose that you view $c$ as a function of, say, the grid size of your finite element discretization of $\Omega$, i.e., $c = c(h)$. You are interested in the limit where $h$ tends to zero. However, I am hopeful, that you can use a large value of $h$ to rapidly locate an interval where $K$ changes sign. This gives you a bracket around a root and you can either apply the bisection method or combine the faster secant method with bisection to obtain a fast and robust method.
Estimating the error on the computed value of $c$ hinges on establishing the existence of an asymptotic error expansion of the form $$c - c(h) = \alpha h^p + \beta h^q + O(h^r), \quad 0 < p < q < r,$$
but that is a topic for another day.
