We apply the Gradient-Newton-Galerkin-Algorithm (GNGA) of Neuberger & Swift to find solutions to a semilinear elliptic Dirichlet problem on the region whose boundary is the Koch snowflake. In a recent paper, we described an accurate and efficient method for generating a basis of eigenfunctions of the Laplacian on this region. In that work, we used the symmetry of the snowflake region to analyze and post-process the basis, rendering it suitable for input to the GNGA. The GNGA uses Newton's method on the eigenfunction expansion coefficients to find solutions to the semilinear problem. This article introduces the bifurcation digraph, an extension of the lattice of isotropy subgroups. For our example, the bifurcation digraph shows the 23 possible symmetry types of solutions to the PDE and the 59 generic symmetry-breaking bifurcations among these symmetry types. Our numerical code uses continuation methods, and follows branches created at symmetry-breaking bifurcations, so the human user does not need to supply initial guesses for Newton's method. Starting from the known trivial solution, the code automatically finds at least one solution with each of the symmetry types that we predict can exist. Such computationally intensive investigations necessitated the writing of automated branch following code, whereby symmetry information was used to reduce the number of computations per GNGA execution and to make intelligent branch following decisions at bifurcation points.2000 Mathematics Subject Classification. 20C35, 35P10, 65N25. 1 2 JOHN M. NEUBERGER, NÁNDOR SIEBEN, AND JAMES W. SWIFTand that the corresponding eigenfunctions {ψ j } j∈N can be chosen to be an orthogonal basis for the Sobolev space H = H 1 0 (Ω) = W 1,2 0 (Ω), and an orthonormal basis for the larger Hilbert space L 2 = L 2 (Ω). The inner products arerespectively (see [1,9,15,17]). Theorem 8.37 and subsequent remarks in [9] imply that the eigenfunctions are in C ∞ (Ω). In [17], properties of the gradients of eigenfunctions near boundary points are explored in light of the lack of regularity of ∂Ω.Using the Gradient-Newton-Galerkin-Algorithm (GNGA, see [26]) we seek approximate solutions u = M j=1 a j ψ j to (1) by applying Newton's method to the eigenfunction expansion coefficients of the gradient ∇J(u) of a nonlinear functional J whose critical points are the desired solutions. The definition of J, the required variational equations, a description of the GNGA, and a brief history of the problem are the subject of Section 2.The GNGA requires as input a basis spanning a sufficiently large but finite dimensional subspace B M = span{ψ 1 , . . . , ψ M }, corresponding to the first M eigenvalues {λ j } M j=1 . As described in [27], a grid G N of N carefully placed points is used to approximate the eigenfunctions. These are the same grid points used for the numerical integrations required by Newton's method. Section 3 briefly describes the process we use for generating the eigenfunctions.Section 4 concerns the effects of symmetry on automated branch following. Th...