Resolve the discrepancy in the 2D first turning point for the cube-domain Bratu equation

Ascertain which numerical estimate is more accurate for the first turning point parameter \lambda^* of the two-dimensional Bratu equation \Delta u + \lambda e^{u} = 0 on the unit square [0,1]^2 with homogeneous Dirichlet boundary conditions—6.808124408 computed via the symmetric finite difference method or 6.808124423 reported by Doedel and Sharifi (2000)—and, thereby, determine the more precise value of \lambda^*.

Background

For the 2D Bratu equation on the unit square, the authors compute a first turning point of 6.808124408 using their symmetric finite difference method, which is extremely close to the earlier value 6.808124423 reported in the literature. They note, however, that their numerical continuation has not fully converged and that the prior paper used a limited number of mesh and collocation points.

Consequently, there is an explicit uncertainty regarding which of the two values is more accurate, motivating a targeted effort to precisely determine \lambda^* in 2D on the cube domain.