Abstract. We seek discrete approximations to solutions u : Ω → R of semilinear elliptic partial differential equations of the form ∆u + fs(u) = 0, where fs is a one-parameter family of nonlinear functions and Ω is a domain in R d . The main achievement of this paper is the approximation of solutions to the PDE on the cube Ω = (0, π) 3 ⊆ R 3 . There are 323 possible isotropy subgroups of functions on the cube, which fall into 99 conjugacy classes. The bifurcations with symmetry in this problem are quite interesting, including many with 3-dimensional critical eigenspaces. Our automated symmetry analysis is necessary with so many isotropy subgroups and bifurcations among them, and it allows our code to follow one branch in each equivalence class that is created at a bifurcation point. Our most complicated result is the complete analysis of a degenerate bifurcation with a 6-dimensional critical eigenspace.This article extends the authors' work in Automated Bifurcation Analysis for Nonlinear Elliptic Partial Difference Equations on Graphs (Int. J. of Bifurcation and Chaos, 2009), wherein they combined symmetry analysis with modified implementations of the gradient Newton-Galerkin algorithm (GNGA, Neuberger and Swift) to automatically generate bifurcation diagrams and solution graphics for small, discrete problems with large symmetry groups. The code described in the current paper is efficiently implemented in parallel, allowing us to investigate a relatively fine-mesh discretization of the cube. We use the methodology and corresponding library presented in our paper An MPI Implementation of a Self-Submitting Parallel Job Queue (Int. J. of Parallel Prog., 2012).