How to treat non-linear term in finite difference solution of $T_{xx}+T_{yy}+aT^2=0$? Can we linearize $T^2$ When solving $\dfrac{\partial^2 T}{\partial x^2}+\dfrac{\partial^2 T}{\partial y^2}+aT^2=0$ by finite difference? 
I solved $\dfrac{\partial^2 T}{\partial x^2}+\dfrac{\partial^2 T}{\partial y^2}=0$ in Matlab using a finite difference explicit scheme. But when there is a source term, I come up with a system of nonlinear algebraic equations and I can't solve it anymore. 
Is there a better method for solving nonlinear equations without linearizing them? 
 A: Just one possible way,
you could use an iterative method for the corresponding parabolic PDE and then treat the nonlinearity as a forcing term.  For example, consider the steady state solution to,
$$ T_{t}(x,y;t) = T_{xx} + T_{yy} + aT^2.
$$ 
Under some assumptions about the sign and possibly magnitude of $a$, you could analytically show that iteratively solving this parabolic PDE (meaning running the iterative method until $T_t \approx 0$) will give you an approximation that converges to the true solution of the original elliptic PDE.  Going on with the assumption that this works, you could discretize using something like Crank Nicolson and let the nonlinearity just be the grid function value at that point squared, i.e. $T^2(x_i,y_j,t_n) = ({\cal T}_{ij}^n)^2,$ where ${\cal T}$ is your grid function or numerical approximation.
A: One family of solutions is $$T(x,y) = -\frac{6 (b^2 + c^2)}{a} \wp(bx + cy + d, 0, e)$$
for arbitrary constants $b,c,d,e$, where $\wp$ is the Weierstrass P function.
A: I would try an iterative approach, solving $\partial_x^2 T_{k+1} + \partial_y^2 T_{k+1} + a T_k^2 = 0$ for $k=0, 1, 2, \ldots$ with $T_0 \equiv 0.$ That is, first solve
$$\dfrac{\partial^2 T_1}{\partial x^2}+\dfrac{\partial^2 T_1}{\partial y^2} = 0,$$
then
$$\dfrac{\partial^2 T_2}{\partial x^2}+\dfrac{\partial^2 T_2}{\partial y^2} + a T_1^2 = 0,$$
$$\dfrac{\partial^2 T_3}{\partial x^2}+\dfrac{\partial^2 T_3}{\partial y^2} + a T_2^2 = 0,$$
and so on until the solution hopefully has converged enough.
A: One simple way to discretize $\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=0$ is
$$\frac{T(x+h,y)+T(x-h,y)-2T(x,y)}{h^2} + \frac{T(x,y+h)+T(x,y-h)-2T(x,y)}{h^2}=0$$
which leads to
$$T(x,y)=\frac{T(x+h,y)+T(x-h,y)+T(x,y+h)+T(x,y-h)}{4}$$
and can be used as an iterative method.
To address nonlinearity, write
$$\frac{T(x+h,y)+T(x-h,y)-2T(x,y)}{h^2} + \frac{T(x,y+h)+T(x,y-h)-2T(x,y)}{h^2}+aT(x,y)^2=0$$
So
$$ah^2T(x,y)^2-4T(x,y)+\left[T(x+h,y)+T(x-h,y)+T(x,y+h)+T(x,y-h)\right]=0$$
This is a quadratic equation on $T(x,y)$. So, you can try iteratively solving this equation.
