  We seek discrete approximations to solutions $u:\Omega \to R$ of semilinear
elliptic partial differential equations of the form $\Delta u + f_s(u) = 0$,
where $f_s$ is a one-parameter family of nonlinear functions and $\Omega$ is a
domain in $R^d$. The main achievement of this paper is the approximation of
solutions to the PDE on the cube $\Omega=(0,\pi)^3 \subseteq 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 {\it 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 {\it An
MPI Implementation of a Self-Submitting Parallel Job Queue} (Int. J. of
Parallel Prog., 2012).
