  We solve by Chebyshev spectral collocation some genuinely nonlinear
Liouville-Bratu-Gelfand type, 1D and a 2D boundary value problems. The problems
are formulated on the square domain $[-1, 1]\times[-1, 1]$ and the boundary
condition attached is a homogeneous Dirichlet one. We pay a particular
attention to the bifurcation branch on which a solution is searched and try to
estimate empirically the attraction basin for each bifurcation variety. The
first eigenvector approximating the corresponding the first eigenfunction of
the linear problem is used as an initial guess in solving the non-linear
algebraic system of Chebyshev collocation to find the "small" solution. For the
same value of the bifurcation parameter we use another initial guess, namely
lowest basis function (1 point approximation), to find the "big" solution. The
Newton-Kantorovich method solves very fast the non-linear algebraic system in
no more than eight iterations. Beyond being exact, the method is numerically
stable, robust and easy to implement. Actually, the MATLAB code essentially
contains three programming lines. It by far surpasses in simplicity and
accuracy various methods used to solve some well-known problems. We end up by
providing some numerical and graphical outcomes in order to underline the
validity and the effectiveness of our method, i.e., norms of Newton updates in
solving the algebraic systems and the decreasing rate of Chebyshev coeffcients
of solution.
