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.