Solving the three-dimensional (3D) Bratu equation is highly challenging due to the presence of multiple and sharp solutions. Research on this equation began in the late 1990s, but there are no satisfactory results to date. To address this issue, we introduce a symmetric finite difference method (SFDM) which embeds the symmetry properties of the solutions into a finite difference method (FDM). This SFDM is primarily used to obtain more accurate solutions and bifurcation diagrams for the 3D Bratu equation. Additionally, we propose modifying the Bratu equation by incorporating a new constraint that facilitates the construction of bifurcation diagrams and simplifies handling the turning points. The proposed method, combined with the use of sparse matrix representation, successfully solves the 3D Bratu equation on grids of up to 3013 points. The results demonstrate that SFDM outperforms all previously employed methods for the 3D Bratu equation. Furthermore, we provide bifurcation diagrams for the 1D, 2D, 4D, and 5D cases, and accurately identify the first turning points in all dimensions. All simulations indicate that the bifurcation diagrams of the Bratu equation on the cube domains closely resemble the well-established behavior on the ball domains described by Joseph and Lundgren [1]. Furthermore, when SFDM is applied to linear stability analysis, it yields the same largest real eigenvalue as the standard FDM despite having fewer equations and variables in the nonlinear system.