\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \usepackage{algorithm} \usepackage{algorithmic} \usepackage{bm} \usepackage[all]{xy} \AtBeginDocument{{\noindent\small Variational and Topological Methods: Theory, Applications, Numerical Simulations, and Open Problems (2012). {\em Electronic Journal of Differential Equations}, Conference 21 (2014), pp. 61--76. ISSN: 1072-6691. http://ejde.math.txstate.edu, http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2014 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \setcounter{page}{61} \title[\hfilneg EJDE-2014/Conf/21 \hfil Symmetry analysis and numerical solutions] {Symmetry analysis and numerical solutions for semilinear elliptic systems} \author[C. Diggans, J. Neuberger, J. Swift \hfil EJDE-2014/Conf/21\hfilneg] {C. Tyler Diggans, John M. Neuberger, James W. Swift} % in alphabetical order \address{C. Tyler Diggans \newline Department of Mathematics, Northern Arizona University, 86005, Flagstaff, AZ 86011, USA} \email{Tyler.Diggans@nau.edu} \address{John M. Neuberger \newline Department of Mathematics, Northern Arizona University, 86005, Flagstaff, AZ 86011, USA} \email{John.Neuberger@nau.edu} \address{Jim W. Swift \newline Department of Mathematics, Northern Arizona University, 86005, Flagstaff, Arizona, USA} \email{Jim.Swift@nau.edu} \thanks{Published February 10, 2014.} \subjclass[2000]{35J15, 65N30} \keywords{Nonlinear elliptic PDE; elliptic systems; Newton's method; \hfill\break\indent GNGA; bifurcation} \begin{abstract} We study a two-parameter family of so-called Hamiltonian systems defined on a region $\Omega$ in $\mathbb{R}^d$ with the bifurcation parameters $\lambda$ and $\mu$ of the form: \begin{gather*} \Delta u + \frac{\partial}{\partial v}H_{\lambda,\mu}(u,v)=0 \quad \text{in } \Omega,\\ \Delta v + \frac{\partial}{\partial u}H_{\lambda,\mu}(u,v)=0, \quad \text{in } \Omega \end{gather*} taking $H_{\lambda,\mu}$ to be a function of two variables satisfying certain conditions. We use numerical methods adapted from \textit{Automated Bifurcation Analysis for Nonlinear Elliptic Partial Difference Equations on Graphs} (Inter. J. Bif. Chaos, 2009) to approximate solution pairs. After providing a symmetry analysis of the solution space of pairs of functions defined on the unit square, we numerically approximate bifurcation surfaces over the two dimensional parameter space. A cusp catastrophe is found on the diagonal in the parameter space where $\lambda=\mu$ and is explained in terms of symmetry breaking bifurcation. Finally, we suggest a more theoretical direction for our future work on this topic. \end{abstract} \maketitle \numberwithin{equation}{section} \allowdisplaybreaks \section{Introduction} Partial Differential Equations (PDE) is an area of mathematics with many applications to physical systems across disciplines. Coupled systems of reaction diffusion equations in particular can be used to model everything from pattern formation of animal coats to chemotaxis in chemical reactions. Our work focuses on numerically approximating steady state solutions of systems of reaction diffusion equations with nonlinear reaction terms that meet certain criteria. In particular, we study two-parameter families of coupled systems of Boundary Value Problems (BVP) with zero Dirichlet boundary conditions of the form: \begin{equation} \label{SYS} \begin{gathered} \Delta u + \frac{\partial}{\partial v}H_{\lambda,\mu}(u,v)=0 \quad \text{in } \Omega, \\ \Delta v + \frac{\partial}{\partial u}H_{\lambda,\mu}(u,v)=0 \quad \text{in } \Omega,\\ u=0,\quad v=0 \quad \text{on }\partial\Omega, \end{gathered} \end{equation} where $\Omega\subseteq \mathbb{R}^d$. Such systems were named Hamiltonian systems in~\cite{Defig}. This is believed to be due to the similarities to Hamilton's equations with two degrees of freedom using the nonconventional representation of kinetic energy, $T=\nabla u\cdot \nabla v$. We take $H_{\lambda,\mu}:\mathbb{R}^2\to \mathbb{R}$ to be a nonlinear function with two real-valued bifurcation parameters $\lambda$ and $\mu$. To apply the methods used in this paper, the function $H_{\lambda,\mu}$ must satisfy three conditions: (1) the point $(0,0)$ is a critical point, (2) the bifurcation parameters $\lambda$ and $\mu$ control the slope of the nonlinearities near the origin, and (3) The quadratic terms of $H_{\lambda,\mu}$ are dependent on either $u$ only or $v$ only. As an example, we will consider the superlinear, subcritical, two-parameter family of systems defined by the choice \begin{equation} H_{\lambda,\mu}(u,v)=\mu\frac{u^2}{2}+\lambda\frac{v^2}{2}+\frac{u^4}{4} +\frac{v^4}{4}.\label{Func} \end{equation} This means that we will be studying the particular nonlinear family of BVP defined by \begin{equation} \label{PDE} \begin{gathered} \Delta u + \lambda v + v^3 =0 \quad \text{in } \Omega, \\ \Delta v + \mu u + u^3 =0 \quad \text{in } \Omega, \\ u=0,\quad v=0 \quad \text{on } \partial\Omega. \end{gathered} \end{equation} From~\cite{Defig}, we know that solution pairs to this Hamiltonian system will be critical points of a related action functional. Letting $H$ represent the Sobelev subspace $H_0^{1,2}(\Omega)$ consisting of $L_2(\Omega)$ functions with one generalized derivative and compact support, we can define the action functional $J_{\lambda,\mu}:H\times H\to \mathbb{R}$ as \begin{equation} \begin{aligned} J_{\lambda,\mu}(u,v) &=\int_{\Omega}\nabla u \cdot \nabla v -H_{\lambda,\mu} \,\mathrm{d}\Omega \\ &=\int_{\Omega}\nabla u \cdot \nabla v -\mu\frac{u^2}{2}-\lambda\frac{v^2}{2} -\frac{u^4}{4}-\frac{v^4}{4} \,\mathrm{d}\Omega. \end{aligned}\label{Functional} \end{equation} We find approximations to solution pairs $U=(u,v)\in H\times H$ using a numerical variational approach similar to that of~\cite{NSS3}. We chose our particular function $H_{\lambda,\mu}$ due to the body of work in existence that deals with related families (see \cite{CCN,NS,NSS3}). In Section~\ref{Prelims}, we discuss the modifications to the numerical methods from~\cite{NSS3} that were needed to apply them to systems of BVP. The contour plots in this paper were made using Mathematica and all other graphics were made using MATLAB. In Section~\ref{Symm}, we present a condensed symmetry digraph relating the possible symmetry types of pairs of functions defined on the unit square $(0,1)^2$. Using the symmetry analysis, we explain the existence of a cusp catastrophe occurring on the $\lambda=\mu$ diagonal in the parameter space in Section~\ref{Surface} . We then show numerical solutions of various symmetry types along with the corresponding portions of the bifurcation diagram in Section~\ref{Results}. Finally, in Section~\ref{Future}, we present a future direction for our research involving an existence proof of a minimal energy sign-changing exactly-once solution pair for the system. \section{Preliminaries} \label{Prelims} We want to find critical points of $J_{\lambda,\mu}$, since they are the solution pairs to \eqref{PDE}. Such points are pairs of functions in $H$, which has an infinite set of basis functions. Each of these basis functions is a continuous function in two variables. Thus, we need to discretize many quantities to represent solutions numerically. This section begins by explaining the methods used to do this. For the convenience of the reader, we then briefly review the Gradient-Newton-Galerkin-Algorithm (GNGA) and its constrained variations from \cite{NSS3}. We focus on the changes made to the algorithms rather than an in-depth explanation of the algorithms themselves. Lastly, we give some simple solutions to\eqref{PDE} by considering a related single BVP. \subsection{Discretization and the Galerkin basis} We take the $M$ lowest energy eigenfunctions of the associated elliptic operator as a basis for a finite subspace of $H$. This type of truncated basis is referred to as a Galerkin basis. The eigenfunctions of the negative Laplacian operator on the region $(0,1)^2$ are the well known products of sine functions and their eigenvalues are the doubly indexed set $\lambda_{m,n}=(m^2+n^2)\pi^2$. We will order these eigenvalues as $\lambda_1=2\pi^2\leq\lambda_2=5\pi^2\leq\lambda_3=5\pi^2\leq\lambda_4=8\pi^2\leq\cdots$ and use a single index for the remainder of this paper. Thus our Galerkin basis for $B_M\subseteq H$ can be defined as \begin{equation*} B_M= \operatorname{span}\big\{\psi_{m,n}=2\sin(m\pi x)\sin(n\pi y): m^2+n^2<(M^*)^2\big\}, \end{equation*} where $M^*$ is chosen to define the number of modes $M$. This ensures that we use a complete set of eigenfunctions of each energy level chosen. For all of the results shown here, unless otherwise specified, we used $M^*=8$, which gives $M=41$ modes. Although this is not a large number of modes, it proved to be sufficient for all the desired results and allowed the speed of the algorithms to be easily performed on a laptop. We divide the region $\Omega=(0,1)^2$ into $N^*$ equal subregions. We define $N=\sqrt{N^*}$ to be the number of subregions in a single dimension, which gives a mesh of $N^2$ evenly spaced grid points on $(0,1)^2$. Using the values of functions at the cell centers, we can numerically represent a function as an $N^2$ column vector. As in \cite{NS}, we choose $N$ large enough to make the numerical integration exact for the products of sine functions required by our algorithm. For all of the plots unless otherwise specified, we in fact use $N=83$ in order to make the contour plots more detailed. We can then approximate solution pairs $U=\left(u,v\right)$ in $H\times H$ as a column vector $\bm{U}\in \mathbb{R}^{2N^2}$. Since our function subspace is now finite dimensional, we can also represent function pairs $U$ as truncated Fourier Sine series. We use a column vector $\bm{c}\in \mathbb{R}^{2M}$, where the entries of $\bm{c}$ are the projections of $\bm{u}$ onto the basis eigenvectors followed by the projections of $\bm{v}$ onto the basis eigenvectors. In particular, we create an $N^2\times M$ matrix $\Psi$ where the $k^{\text{th}}$ column of $\Psi$ is the vector approximation $\bm{\psi_k}\in \mathbb{R}^{N^2}$ of the eigenfunction $\psi_k$ corresponding to the eigenvalue $\lambda_k$. We then have the relationship \[ \bm{U}=\begin{bmatrix} \bm{u}\\ \bm{v} \end{bmatrix}=[\Psi| \Psi] \bm{c}. \] \subsection{Derivatives of the action functional} As known, we can find solutions to system~\eqref{PDE} by finding critical points of the related functional \eqref{Functional}. We consider a pair of functions $U$ to be a critical point of $J_{\lambda,\mu}$ if the directional derivative of $J_{\lambda,\mu}$ at that point is zero in all of the Galerkin basis directions. Using the limit definition of the derivative, the directional derivative in a general direction $(w,z)$ is given by \begin{align*} J_{\lambda,\mu}'(u,v)(w,z) &=\lim_{t\to 0}\frac{1}{t}\Big[\int_\Omega\nabla (u+tw) \cdot \nabla (v+tz) - \mu \frac{(u+tw)^2}{2} -\lambda \frac{(v+tz)^2}{2} \\ &\quad -\frac{(u+tw)^4}{4} -\frac{(v+tz)^4}{4} \,\mathrm{d}\Omega\\ &\quad - \int_\Omega\nabla u \cdot \nabla v - \mu \frac{u^2}{2} -\lambda \frac{v^2}{2} -\frac{u^4}{4} -\frac{v^4}{4} \phantom{x} \mathrm{d}\Omega\Big] \\ &=\int_\Omega\nabla u \cdot \nabla z +\nabla w \cdot \nabla v - \mu u w -\lambda v z -u^3w -v^3z \, \mathrm{d}\Omega. \end{align*} Since our Galerkin basis is orthonormal, it is sufficient to check that the directional derivatives in the $(\psi_k,0)$ and $(0,\psi_k)$ directions are zero for each $k\in\{1,2,\dots,M\}$. One can perform calculus such as this for all of $J_{\lambda,\mu}$'s various partial derivatives for general nonlinearities $H_{\lambda,\mu}$. For our choice of $H_{\lambda,\mu}$, the derivative in the $(\psi_k,0)$ direction simplifies to \begin{equation} J_{\lambda,\mu}'(u,v)(\psi_k,0) =\int_\Omega{\nabla \psi_k \cdot \nabla v -\mu u \psi_k - u^3 \psi_k \phantom{x} \,\mathrm{d}\Omega}. \label{DirD} \end{equation} The derivative in the $(0,\psi_k)$ direction is defined similarly. We will also need the second directional derivatives of $J_{\lambda,\mu}$ in various basis directions in order to implement Newton's method on a vector function of the first derivatives of $J_{\lambda,\mu}$. The general result in the $((\psi_i,0),(\psi_j,0))$ direction is given by \[ J_{\lambda,\mu}''(u,v)(\psi_i,0)(\psi_j,0) =\int_\Omega{-\mu\psi_i\psi_j-3 u^2 \psi_i \psi_j \mathrm{d}\Omega}. \] Derivatives in the $((0,\psi_i),(0,\psi_j))$ directions are defined similarly. The mixed directional derivatives in the $((\psi_i,0),(0,\psi_j))$ direction or the $((0,\psi_i),(\psi_j,0))$ direction both simplify to \[ \int_\Omega \nabla \psi_i \cdot \nabla \psi_j \mathrm{d}\Omega =\lambda_{i}\delta_{i,j}, \] where $\delta_{i,j}$ is the Kronecker delta function. The next subsection uses the derivatives defined above to describe the details of the numerical variational approach adapted from \cite{NSS3} and used in this paper. \subsection{Numerical methods} \label{Numerical} We will be performing a $2M$ dimensional Newton's method to find zeros of a vector function $g_{\lambda,\mu}:\mathbb{R}^{2M}\to \mathbb{R}^{2M}$ that approximates the directional derivatives of $J_{\lambda,\mu}$ in the Galerkin basis directions. We use Green's second theorem and the orthogonality of the basis functions to simplify the leading integration terms in \eqref{DirD} to \[ \int_\Omega \nabla \psi_i \cdot \nabla v -\mu u \psi_i \mathrm{d}\Omega = -\lambda_i c_{M+i} -\mu c_i. \] We then use a Riemann sum to approximate the remaining term of the integral by \[ \int_\Omega u^3 \psi_i \mathrm{d}\Omega \approx \frac{\bm{\psi_i}^T \bm{u}^3}{N^2}, \] where $\bm{u}^3$ represents $\bm{u}$ cubed component-wise. Thus, we are finding critical points of an approximated gradient vector function $g_{\lambda,\mu}$ of $J_{\lambda,\mu}$, where the $i^{\textrm{th}}$ component of $g_{\lambda,\mu}$ is given by \[ g_{\lambda,\mu}(\bm{c})_i =\begin{cases} -\lambda_i c_{i+M}-\mu c_i -\frac{\psi_i^T \bm{u}^3}{N^2} & \text{if } i\leq M \\ -\lambda_{i-M} c_{i-M}-\lambda c_i -\frac{\psi_{i-M}^T \bm{v}^3}{N^2} & \text{if } i>M. \end{cases} \] To implement Newton's method to find a zero of $g_{\lambda,\mu}$, we require the Jacobian of the vector function $g_{\lambda,\mu}$. For this we use the approximated Hessian matrix of $J_{\lambda,\mu}$. Using numerical integration as above, the elements of $h_{\lambda,\mu}$ are calculated as \[ h_{\lambda,\mu}(\bm{c}) =\left[\begin{array}{c|c} \\ \Big(-\mu \delta_{ij} - 3 \frac{(\bm{\psi}_i\cdot\bm{\psi}_j)^T \bm{u}^2}{N^2}\Big)_{ij} & \Lambda \\ \\ \hline \\ \Lambda & \Big(-\lambda \delta_{ij} - 3\frac{(\bm{\psi}_i\cdot\bm{\psi}_j)^T \bm{v}^2}{N^2}\Big)_{ij}\\ \\ \end{array}\right], \] where $\Lambda$ is the diagonal matrix of the ordered eigenvalues of the basis functions of $B_M$. Using $g_{\lambda,\mu}$ and $h_{\lambda,\mu}$ as described above, the main algorithm for GNGA is the same as in \cite{NSS3}. For the convenience of the reader we include it in Algorithm \ref{GNGA}: \begin{algorithm} \caption{(GNGA) Gradient-Newton-Galerkin Algorithm.} \label{GNGA} \begin{algorithmic}[1] \REQUIRE An initial guess $\bm{c}_{\rm guess}$, tol \ENSURE $\bm{c}$ such that $g_{\lambda,\mu}(\bm{c})=0$ \STATE Set $\bm{c}=\bm{c}_{\rm guess}$\\ \STATE Compute $g:=g_{\lambda,\mu}(\bm{c})$\\ \WHILE{$|g|>\textrm{tol}$} \STATE Compute $h:=h_{\lambda,\mu}(\bm{c})$\\ \STATE Solve $g=h\chi$ for search direction $\chi$ \\ \STATE Set $\bm{c}=\bm{c}-\chi$\\ \STATE Update $g:=g_{\lambda,\mu}(\bm{c})$\\ \ENDWHILE \end{algorithmic} \end{algorithm} This algorithm only differs from the one used in \cite{NSS3} by the definitions of $g_{\lambda,\mu}$ and $h_{\lambda,\mu}$. Since the methods in~\cite{NSS3} were designed to follow bifurcation branches, we restrict our searches in the two-dimensional parameter space to one-dimensional diagonal lines of the form $\lambda=\mu+s$, with $s$ fixed in $\mathbb{R}$. In this way, we have only one bifurcation parameter when following a particular branch. As in~\cite{NSS3}, we use two different constrained GNGA methods to find and follow new branches for a fixed $s$ value. The appended vector $\tilde{\bm{c}}=(\bm{c},\lambda)$ is considered a solution if $g_{\lambda,\mu}(\bm{c})=0$ for the parameter value $\lambda$. The constrained GNGA methods allow the bifurcation parameter $\lambda$ to vary as an additional variable. This brings the need for an additional equation in order to have a determined system. In each constrained variation of GNGA we define a constraint equation $\kappa\left(\tilde{\bm{c}}\right)=0$. We define an appended gradient vector $\widetilde{g}_{\lambda,\mu}$ and Hessian matrix $\widetilde{h}_{\lambda,\mu}$ for a given $\kappa$ as \[ \widetilde{g}_{\lambda,\mu}(\bm{c}) =\left[\begin{array}{c} g_{\lambda,\mu}(\bm{c})\\ \\ \hline 0 \end{array}\right] \quad \text{ and }\quad \widetilde{h}_{\lambda,\mu}(\bm{c}) =\left[ \begin{array}{c|c} h_{\lambda,\mu}(\bm{c})&\frac{\partial g_{\lambda,\mu}}{\partial \lambda}\\ &\\ \hline (\nabla_{\bm{c}}\phantom{x} \kappa)^T& \frac{\partial\kappa}{\partial \lambda} \end{array}\right]. \] Appending the zero to $g_{\lambda,\mu}$ ensures $\tilde{\bm{c}}$ satisfies the constraint equation $\kappa(\tilde{\bm{c}})=0$. To implement the constrained variations of GNGA, we use Algorithm~\ref{GNGA} with the substitution of $\widetilde{g}_{\lambda,\mu}$ and $\widetilde{h}_{\lambda,\mu}$ for $g_{\lambda,\mu}$ and $h_{\lambda,\mu}$ respectively. We now describe the numerical approach in the order in which our code implements each stage. We start by following the known branch of trivial solutions. When at a bifurcation point, we find and follow new branches. The tangent-GNGA (tGNGA) method is used to follow a symmetry invariant branch, given an old and a current solution denoted $\tilde{\bm{c}}_{\textrm{old}}$ and $\tilde{\bm{c}}_{\textrm{cur}}$. We begin by computing an approximated tangent vector as \[ \bm{T}=(\tilde{\bm{c}}_{\textrm{cur}}-\tilde{\bm{c}}_{\textrm{old}})/ \|\tilde{\bm{c}}_{\textrm{cur}}-\tilde{\bm{c}}_{\textrm{old}}\|\in \mathbb{R}^{2M+1}. \] We then choose an initial guess by setting $\tilde{\bm{c}}_{\rm guess}=\tilde{\bm{c}}_{\textrm{cur}}+b \bm{T}$ where $b$ is the branch-following speed. Here we used a less dynamic version of the algorithm than the one given in~\cite{NSS3}, in which our step size $b$ is fixed at $b=0.5$. When more detail is needed, we change the step size manually. We define the tangent constraint to be $\kappa_t(\tilde{\bm{c}}):=(\tilde{\bm{c}}-\tilde{\bm{c}}_{\rm guess})\cdot \bm{T}$. Hence the result $\tilde{\bm{c}}=\textbf{tGNGA}(\tilde{\bm{c}}_{guess},\bm{T})$ satisfying $\kappa_t(\tilde{\bm{c}})=0$ lies in the hyperplane passing through our guess and perpendicular to $\bm{T}$. See~\cite{NSS3} for more details on this method. From the implicit function theorem, we know that bifurcation points can only occur when the Hessian matrix has a zero eigenvalue. Thus, while following a symmetry invariant branch, we monitor the Morse index (MI) of the approximated Hessian matrix $h_{\lambda,\mu}$. When a change in Morse index is observed between two consecutive solutions, say from a MI of $k$ at $\tilde{\bm{c}}_{\textrm{old}}$ to a MI of $k+\delta$ at $\tilde{\bm{c}}_{\textrm{cur}}$, we implement a vector secant method to find the bifurcation point(s) between them. For the particular $k$ value we first define a function $\beta:\mathbb{R}^{2M+1}\to\mathbb{R}$ that maps solution vectors $\tilde{\bm{c}}$ to the $k^\textrm{th}$ eigenvalue of $h_{\lambda,\mu}(\bm{c})$. Then given the two solutions $\tilde{\bm{c}}_0$ and $\tilde{\bm{c}}_1$ where the MI changes and the fixed, approximate tangent vector $T=(\tilde{\bm{c}}_1-\tilde{\bm{c}}_0)/\|\tilde{\bm{c}}_1-\tilde{\bm{c}}_0\|$, we iterate the following until $\beta(\tilde{\bm{c}}_i)$ is less than a given tolerance: \begin{itemize} \item Set $\tilde{\bm{c}}_{\rm guess}=\tilde{\bm{c}}_i -\frac{(\tilde{\bm{c}}_i-\tilde{\bm{c}}_{i-1}) \beta(\tilde{\bm{c}}_i)}{(\beta(\tilde{\bm{c}}_i) -\beta(\tilde{\bm{c}}_{i-1}))}$ \item Set $\tilde{\bm{c}}_{i+1}=\textrm{\textbf{tGNGA}} (\tilde{\bm{c}}_{\rm guess},\bm{T})$ \item Set $i=i+1$. \end{itemize} Once we find the bifurcation point $\tilde{\bm{c}}^*$ using the vector secant method, we then use the cylinder-GNGA (cGNGA) method adapted from \cite{NSS3} to find a solution with a nonzero projection onto a subspace $E$ of the critical eigenspace of $h_{\lambda,\mu}(\tilde{\bm{c}}^*)$. To do this we constrain our search to solutions that lie on the cylinder $C=\{\tilde{\bm{c}} \in \mathbb{R}^{2M+1} : \|P_E(\tilde{\bm{c}}-\tilde{\bm{c}}^*)\|=\epsilon\}$, where $P_E$ is the orthogonal projection onto $E$ and $\epsilon$ is a small fixed parameter. This requirement leads us to the constraint $\kappa_c(\tilde{\bm{c}})=\frac{1}{2}(\|P_E(\tilde{\bm{c}}-\tilde{\bm{c}}^*)\|^2-\epsilon^2)=0$. We use an initial guess of $\tilde{\bm{c}}_{\rm guess}=\tilde{\bm{c}}^*+\epsilon\hat{e}$, where $\hat{e}$ is a randomly chosen unit vector in the space $E$. We then use the bifurcation point $\tilde{\bm{c}}^*$ and the result $\tilde{\bm{c}}_\textrm{cur}=\textbf{cGNGA}(\tilde{\bm{c}}^*)$ to implement $tGNGA$ to follow this new branch. In this way, we can theoretically find all solutions connected to the trivial bifurcation surface. The next subsection gives our initial solutions to Equation~\eqref{PDE} by considering a special case. \subsection{A related single parameter family of PDE} Our first solutions to \eqref{PDE} come from the restriction of setting $\lambda=\mu$ and $u=v$. In this case, our system becomes two copies of a well known family of single BVP given by \begin{equation} \label{Single} \begin{gathered} \Delta u + \lambda u + u^3= 0 \quad\text{in } (0,1)^2 \\ u=0,\quad v=0 \quad \text{on }\partial (0,1)^2 . \end{gathered} \end{equation} This family was studied in \cite{NS} in detail. Thus, we know that the trivial function $u\equiv 0$ is a solution to \eqref{Single} for all values of $\lambda$, and that there are primary branches bifurcating from the trivial branch when $\lambda$ equals the eigenvalues of the elliptic operator (i.e. $\lambda=\lambda_{m,n}=(m^2+n^2)\pi^2$ for $m,n\in\mathbb{Z}^+$). Note that there are no bifurcations from the trivial branch for negative parameter values in this restrictive case. We will find much richer bifurcation behavior in the system including bifurcations for negative parameter values, but the solutions given here give us an intuition of what types of solution pairs to expect in the system. \section{Symmetry of the solution space} \label{Symm} We present a symmetry analysis of the solution space consisting of pairs of functions defined on the unit square. The study of any dynamical system that can be represented as a coupled system of PDE defined on a pair of unit squares will benefit from the symmetry analysis in this section. In symmetry breaking bifurcations, we generally see solution branches of more symmetry bifurcating to solution branches of less symmetry. Thus, knowing the full symmetry digraph of a solution space will allow us to predict and understand the connections between parent and child branches in the bifurcation diagram. In addition, we can assure whether solutions of all possible symmetry types have been found for a given system. The symmetry of a pair of functions defined on the unit square is isomorphic to the group $D_4\times\mathbb{Z}_2\times \mathbb{Z}_2$. Figure \ref{SymmPic} shows a pair of unit squares in gray along with the generators chosen to represent the symmetry group of $D_4\times\mathbb{Z}_2$. Here $\rho$ is the $90\deg$ rotation of both functions, $\tau$ is the reflection of both functions about a particular diagonal, and $\sigma$ is the symmetry between the two functions. The second copy of $\mathbb{Z}_2$ takes into account that a function may be positive or negative at a given point. This is represented by the generator $\langle -1\rangle$ \begin{figure}[ht] \begin{center} \includegraphics[width=0.4\textwidth]{fig1} % SymmPic \end{center} \caption{A pair of unit squares whose symmetry is isomorphic to $D_4\times\mathbb{Z}_2$ and can be generated by $\langle \rho, \tau, \sigma\rangle$} \label{SymmPic} \end{figure} From \cite{NSS3}, we know that $D_4$ symmetry is isomorphic to the symmetry of the graph shown in Figure~\ref{Graph} (a). Thus when studying a pair of functions on the unit square, the symmetry of the pair can be modeled by two copies of this graph, where the two vertices of each corner of the squares are identified. This graph can be represented in three dimensions as a polygon with two octagonal faces, four hexagonal faces and eight triangles as shown in Figure~\ref{Graph} (b). \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}} \includegraphics[width=0.3\textwidth]{fig2a} \quad % D4Graph &\includegraphics[width=0.35\textwidth]{fig2b} % perspectivegraph \\ (a) & (b) \end{tabular} \end{center} \caption{(a) A graph with the same symmetry as $D_4$. (b) A graph with the same symmetry as $D_4 \times \mathbb{Z}_2$} \label{Graph} \end{figure} The group $D_4\times \mathbb{Z}_2\times \mathbb{Z}_2=\langle \rho, \tau, \sigma, -1\rangle$ has $92$ different symmetry types. We use the GAP program to condense these to a more manageable $20$ condensation classes of symmetry types. See~\cite{GAP} for more information on this program. To illustrate the need for this condensation, the output symmetry digraphs from the GAP program are included in Figure~\ref{Digraphs}. \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}} \includegraphics[width=0.45\textwidth]{fig3a} % digraph92 &\includegraphics[width=0.45\textwidth]{fig3b} % digraph20 \\ (a) & (b) \end{tabular} \end{center} \caption{The two output digraphs from the GAP program \cite{GAP}. The full symmetry digraph of $D_4\times \mathbb{Z}_2\times \mathbb{Z}_2$ is given in (a), while the condensed symmetry digraph of $D_4\times \mathbb{Z}_2\times \mathbb{Z}_2$ is given in (b)} \label{Digraphs} \end{figure} GAP identifies two symmetry types $A$ and $B$ as related, if there is an automorphism of $D_4 \times \mathbb{Z}_2\times \mathbb{Z}_2$ mapping a representative of $A$ to a representative of $B$. The automorphisms of our group can be described as a $45$ degree rotation of the graph in Figure~\ref{Graph}(b), and replacing any generator of a subgroup with its negation. All symmetry types that are related by one of these automorphisms make up what is referred to as a condensation class. The symmetry digraph of condensation classes is shown with more detail in Figure~\ref{CONDDIGRAPH}. We will use the equivalence class notation $[\langle \rho, \tau, \sigma\rangle]$ to represent the condensation class for which $\langle \rho,\tau, \sigma\rangle$ is a member. For easy reference we will let $\mathcal{C}_i$ denote the $i^{\textrm{th}}$ condensation class. For example the condensation class $\mathcal{C}_1=[\langle \rho, \tau, \sigma\rangle]$ consists of the eight symmetry types of the form $\langle \pm\rho, \pm\tau, \pm\sigma\rangle$. %\begin{landscape} {\tiny %\small \begin{figure}[ht] \xymatrixrowsep{.01in} %\xymatrixrowsep{.02in} \xymatrixcolsep{.01in} \centerline{ \xymatrix{ &&&&&&\mathcal{C}_0&&&&&&&\\ &&&&&&[\langle\rho,\tau,\sigma,-1\rangle]\ar@{->}[ddddlll] \ar@{->}[ddddddddddrrrrr]^{D_4}&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&\mathcal{C}_1&&&&&&&&&&\\ &&&[\langle\rho,\tau,\sigma,\rangle]\ar@{->}[dddddll]\ar@{->}[dddddl] \ar@{->}[ddddd]\ar@{->}[dddddr]\ar@{->}[dddddddddddrrrrr]^(.65){D_4}&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &\mathcal{C}_2&\mathcal{C}_6&\mathcal{C}_5&\mathcal{C}_4&&&&&&&\mathcal{C}_3&&\\ &[\langle\rho,\sigma\rangle]\ar@{->}[dddddl]\ar@{->}[ddddd] \ar@{->}[dddddddddddr]^(.27){\mathbb{Z}_4}&[\langle\rho^2,\rho\tau, \sigma\rangle]\ar@{->}[ddddd]\ar@{->}[dddddll]\ar@{->}[dddddrrrrrr] \ar@{->}[dddddr]&[\langle\rho,\tau\rangle]\ar@{->}[dddddll] \ar@{->}[ddddd]\ar@{->}[dddddl]\ar@{->}[dddddddddddrr]^(.52){D_4}&[\langle\rho^2, \rho\tau,\rho\sigma\rangle]\ar@{->}[dddddlll]\ar@{->}[dddddll] \ar@{->}[dddddl]\ar@{->}[dddddddddddr]^{D_4} \ar@{->}[dddddddddddrrrr]^(.53){D_4}&&&&&&&[\langle-\rho^2,\tau,\sigma \rangle]\ar@{->}[dddddlll]\ar@{->}[dddddl]\ar@{->}[ddddd]\ar@{->}[dddddr]&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ \mathcal{C}_{12}&\mathcal{C}_{13}&\mathcal{C}_9&\mathcal{C}_{11}&&&&& \mathcal{C}_{14}&&\mathcal{C}_7&\mathcal{C}_{10}&\mathcal{C}_8&\\ [\langle\rho^2,\sigma\rangle]\ar@{->}[dddddr]\ar@{->}[dddddrr]&[\langle\rho \rangle]\ar@{->}[ddddd]\ar@{->}[ddddddddddrrrrr]^(.55){\mathbb{Z}_4}&[\langle \rho^2,\rho\tau\rangle]\ar@{->}[ddddd]\ar@{->}[dddddrrr]&[\langle\rho^2, \rho\tau\sigma\rangle]\ar@{->}[dddddll]\ar@{->}[dddddrr]\ar@{->}[dddddrrrrr]&&&&&[ \langle\rho\tau,\sigma\rangle]\ar@{->}[ddddd]\ar@{->}[dddddlll] \ar@{->}[dddddllllll]&&[\langle-\rho^2,\sigma\rangle]\ar@{->}[ddddd] \ar@{->}[dddddllllllll]&[\langle-\rho^2,\tau\rangle]\ar@{->}[dddddllllll] \ar@{->}[dddddl]&[\langle-\rho^2,\rho\tau\sigma\rangle]\ar@{->}[dddddllll] \ar@{->}[dddddll]&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &\mathcal{C}_{18}&\mathcal{C}_{16}&&&\mathcal{C}_{19}&&&\mathcal{C}_{17}&& \mathcal{C}_{15}&&&\\ &[\langle\rho^2\rangle]\ar@{->}[ddddrrrrr]&[\langle\sigma\rangle] \ar@{->}[ddddrrrr]&&&[\langle\rho\tau\rangle]\ar@{->}[ddddr]&&&[ \langle\rho\tau\sigma\rangle]\ar@{->}[ddddll]&&[\langle-\rho^2\rangle] \ar@{->}[ddddllll]&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&&&&&&&&\\ &&&&&&\mathcal{C}_{20}&&&&&&&\\ &&&&&&[\langle 1\rangle]&&&&&&&\\ } } \caption{The Hasse diagram of condensation classes of symmetry types of the solution space of Equation~\eqref{PDE}. We use the usual class notation $[\langle \rho, \tau, \sigma\rangle]$ to denote the set of all symmetry types that are related to $\langle \rho, \tau, \sigma\rangle$ through an automorphism of the group. All directed edges represent a loss of $\mathbb{Z}_2$ symmetry unless otherwise specified} \label{CONDDIGRAPH} %fig4 \end{figure}} %\end{landscape} \section{Bifurcation surfaces and the cusp catastrophe} \label{Surface} To create bifurcation surfaces over the two dimensional parameter space, we begin by observing that the trivial pair $U=(\bm{0},\bm{0})$ is a solution to Equation~\eqref{PDE} for all choices of the parameters $\lambda$ and $\mu$. Hence, the trivial plane is our first bifurcation surface. As discussed in Section~\ref{Prelims}, we know that bifurcations (and turning points) are found when the approximated Hessian matrix $h_{\lambda,\mu}$ is not invertible. For any set values of $\lambda$ and $\mu$ in the trivial plane, the integration terms in $h_{\lambda,\mu}$ vanish and the Hessian becomes the block diagonal matrix \[ h_{\lambda,\mu}(\bm{0})=\left[ \begin{array}{c|c} \text{diag}(-\mu) & \lambda\\ \\ \hline \\ \lambda & \text{diag}(-\lambda)\\ \end{array} \right]. \] This matrix is non invertible exactly when the following $2\times 2$ matrix has a determinant of zero, for some $k\in\left\{1,2,\dots,M\right\}$: \[ \begin{bmatrix} -\mu &\lambda_k \\ \lambda_k& -\lambda \end{bmatrix} \] This occurs when $\lambda\mu=\lambda_k^2$. Thus, for each $k$, we get two such hyperbolas in the trivial plane: one in which both $\lambda$ and $\mu$ are positive and one in which both are negative. These critical hyperbolas are the primary bifurcation curves. Note that bifurcation curves in the two parameter case are analogous to bifurcation points in the one parameter case. In theory, we will have infinitely many positive and negative primary bifurcation surfaces. Numerically, we find $2M$ primary bifurcation surfaces, half of which bifurcate in the positive parameter space and half in the negative parameter space. Figure $5$ shows the critical hyperbolas of the form $\lambda\mu=\lambda_k^2$ for $k\leq 8$. \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig5} %Prim \end{center} \caption{The 10 critical hyperbolas of the form $\lambda\mu=\lambda_k^2$ for $k\leq8$. These are the first five bifurcation curves in the positive and negative parameter spaces. From these curves, $16$ primary bifurcation surfaces originate} \label{Primary} \end{figure} Note that the other two quadrants of the parameter space do not have any bifurcation since the product $\lambda\mu<0$ and $\lambda_k^2>0$. We know that for $\lambda=\mu>0$, we have the primary branches consisting of identical pairs of solutions to the single PDE~\eqref{Single}. In the negative parameter space, the restriction of setting $u=-v$ gives the system \begin{gather*} \Delta u -\mu u - u^3 =0\\ -\Delta u +\mu u + u^3 =0, \end{gather*} which is also two copies of a single PDE. Note that here $\mu<0$, allowing the linear term to dominate the concavity of the solution for small $u$. Thus, $u$ can be positive and concave down satisfying the zero Dirichlet boundary conditions. In the more general case when $\lambda\neq\mu$ we find that $u$ and $v$ are not as simply related. On the line $\lambda=\mu+s$, we find that the primary bifurcation points are found when $\lambda(\lambda-s)=\lambda_k^2$. The quadratic formula gives the $k^{\text{th}}$ pair of bifurcation points occurring when \[ \lambda=\frac{s\pm \sqrt{s^2+4 \lambda_k^2}}{2}. \] It is important to notice that the full richness of symmetry in our system will only be found on the $\lambda=\mu$ diagonal since no other solution pairs will be invariant under $\langle\pm\sigma\rangle$. Thus, when $s\neq 0$ the symmetry breaking bifurcations will not follow the digraph given in Figure~\ref{CONDDIGRAPH}. The symmetry behavior of these types of solutions will follow that given in~\cite{NS}, with the difference of two unrelated functions on the unit square instead of the single function discussed therein. This fundamental difference in symmetry was first encountered when attempting to create the primary bifurcation surfaces. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig6} %Cusp_60deg \end{center} \caption{The first positive primary bifurcation surface using only diagonals where $s\neq0$ and beginning at the critical hyperbola for $\lambda_1=2\pi^2$. We show the norm of the solution vector $\bm{U}$ plotted against the $\lambda$, $\mu$ plane. For reference, the first primary bifurcation curve and the first secondary bifurcation curve on the $\lambda=\mu$ diagonal are included as dotted and dashed-dotted lines respectively. Note the cusp catastrophe that occurs on the $\lambda=\mu$ diagonal. The surface was made using $S=187.5$ and increments of $15$ in order to avoid including the $s=0$ diagonal. Note that this primary surface follows the secondary bifurcation on the $\lambda=\mu$ line, but does not bifurcate for other diagonals} \label{Surf} \end{figure} To approximate a bifurcation surface, we create polygons using linear interpolation of solutions from diagonal lines for varying values of $s$ ranging from some $-S$ to $S$. Figure~\ref{Surf} shows the primary bifurcation surface created for the first positive critical hyperbola. The primary bifurcation curve for the $s=0$ diagonal is included as a dotted line and the first secondary bifurcation curve for $s=0$ is shown as a dashed-dotted line for reference. Notice that the surface appears to follow just above the first secondary bifurcation curve on the $\lambda=\mu$ diagonal. This first secondary bifurcation on the $s=0$ line corresponds to a loss of invariance under the transformation $\langle \sigma\rangle$, thus any other diagonal such that $s\neq0$ will not have this bifurcation since its primary branch will not be invariant under that transformation either. This phenomena is known as a cusp catastrophe. The pitchfork bifurcation present on the $\lambda=\mu$ diagonal becomes a disconnected cusp bifurcation for all other diagonal lines in the parameter space. Figure~\ref{Cusp}(a) shows the usual pitchfork bifurcation found when $s=0$. A contour plot of the solution pair corresponding to the black diamond in (a) is given in Figure~\ref{Cusp}(b). Figure~\ref{Cusp}(c) shows the cusp bifurcation for the diagonal when $s=5$. Note that the two pieces of the diagram are disconnected. We were able to find the curve bifurcating from infinity by using the $s=0$ solution pair as an initial guess in the GNGA method with $\lambda=\mu+5$. Figure~\ref{Cusp}(d) is a contour plot of the solution pair corresponding to the black diamond in (c). \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}c@{}c@{}} \includegraphics[width=0.32\textwidth]{fig7a} % Cusp1.png &\includegraphics[width=0.12\textwidth]{fig7b} % Cusp_Soln_1.png &\includegraphics[width=0.32\textwidth]{fig7c} % Cusp5.png &\includegraphics[width=0.12\textwidth]{fig7d} % Cusp_Soln_5.png \\ (a) &(b) &(c) &(d) \end{tabular} \end{center} \caption{(a) The bifurcation diagram of the first primary curve on the $\lambda=\mu$ diagonal including the secondary bifurcation corresponding to a loss of $\langle \sigma\rangle$ symmetry. Here and in subsequent figures, open circles will be used to denote bifurcation points. (b) The contour plot of the solution pair corresponding to the diamond in (a). (c) The bifurcation diagram of the first primary curve on the $\lambda=\mu+5$ diagonal lacking the secondary bifurcation seen in (a), together with the disconnected portion that bifurcates from infinity. (d) The contour plot of the solution pair corresponding to the diamond in (c)} \label{Cusp} \end{figure} In Figure~\ref{CUSPERS}, we also include an alternative plotting scheme for the vertical axis to display the cusp bifurcation more clearly. Contour plots of the solution pairs marked by the black diamonds from top to bottom are also included. \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}c@{}c@{}} \includegraphics[width=0.42\textwidth]{fig8a} % CUSPERS.png &\includegraphics[width=0.15\textwidth]{fig8b} % CuspSoln1.png &\includegraphics[width=0.15\textwidth]{fig8c} % CuspSoln2.png &\includegraphics[width=0.15\textwidth]{fig8d} % CuspSoln3.png \\ (a) & (b) & (c) & (d) \end{tabular} \end{center} \caption{(a) A plot of $c_5-c_{M+5}$ versus the bifurcation parameter $\lambda$. This shows the cusp bifurcation behavior more clearly. The contours in (b), (c), and (d) show the solutions marked by the black diamonds in (a) from top to bottom} \label{CUSPERS} \end{figure} We now connect our numerical results for the diagonal where $\lambda=\mu$ with the symmetry analysis in Section~\ref{Symm}. \section{Symmetry results for the diagonal where $\lambda=\mu$}\label{Results} The trivial solution pair is invariant under all possible symmetry transformations including $\langle -1 \rangle$. Thus, it represents the very top vertex of the symmetry digraph in Figure~\ref{CONDDIGRAPH}, labeled as $\mathcal{C}_0$. Non degenerate bifurcations from the trivial branch lead to solution pairs that are invariant under symmetry types found in the condensation classes $\mathcal{C}_1$ and $\mathcal{C}_3$ only. We also have found some degenerate bifurcations from the trivial branch, such as when $\lambda=\mu=10\pi^2$ that lead to solution pairs invariant under symmetry types found in $\mathcal{C}_6$. It is important to note that on this diagonal, all primary branches are invariant under either $\langle\sigma\rangle$ or $\langle-\sigma\rangle$. Figures \ref{Prims} and \ref{NegPrim} show bifurcation diagrams and the accompanying contour plots and symmetry types of solution pairs for the first four primary bifurcation branches where $\lambda=\mu$ for positive parameter values and negative parameter values, respectively. \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}c@{}c@{}c@{}} \includegraphics[width=0.37\textwidth]{fig9a} % Prims.png &\includegraphics[width=0.14\textwidth]{fig9b} % Prim_1_1.png &\includegraphics[width=0.14\textwidth]{fig9c} % Prim_1_2_2.png &\includegraphics[width=0.14\textwidth]{fig9d} % Prim_1_2_1.png &\includegraphics[width=0.14\textwidth]{fig9e} % Prim_2_2.png \\ & $\langle \rho,\tau,\sigma\rangle$ & $\langle -\rho^2,\tau,\sigma\rangle$ & $\langle -\rho^2,-\rho\tau,\sigma\rangle$ & $\langle -\rho,\tau,\sigma\rangle$ \\ (a) & (b) & (c) & (d) & (e) \end{tabular} \end{center} \caption{(a) A plot of the norm of the solution vector $\bm{U}$ versus $\lambda$ for the first four positive primary curves on the $\lambda=\mu$ diagonal. (b), (c), (d), and (e) give the contour plots of the solution pairs marked by the black diamonds in (a) from bottom to top.} \label{Prims} \end{figure} \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}c@{}c@{}c@{}} \includegraphics[width=0.37\textwidth]{fig10a} % NPrims.png &\includegraphics[width=0.14\textwidth]{fig10b} % Prim_n1_n1.png &\includegraphics[width=0.14\textwidth]{fig10c} % Prim_n1_n2_2.png &\includegraphics[width=0.14\textwidth]{fig10d} % Prim_n1_n2_1.png &\includegraphics[width=0.14\textwidth]{fig10e} % Prim_n2_n2.png \\ & $\langle \rho,\tau,-\sigma\rangle$ & $\langle -\rho^2,-\tau,-\sigma\rangle$ & $\langle -\rho^2,\rho\tau,-\sigma\rangle$ & $\langle -\rho,\tau,-\sigma\rangle$ \\ (a) & (b)& (c) & (d) & (e) \end{tabular} \end{center} \caption{(a) A plot of the norm of the solution vector $\bm{U}$ versus $\lambda$ for the first four negative primary curves on the $\lambda=\mu$ diagonal. (b), (c), (d), and (e) give the contour plots of the solution pairs marked by the black diamonds in (a) in the order (d), (e), (c), and (b) from top to bottom} \label{NegPrim} \end{figure} Notice that the solution pairs along the branches bifurcating from the first and third bifurcation points are members of $\mathcal{C}_1$ whereas those bifurcating from the second are members of $\mathcal{C}_3$. Solution pairs on the fifth and sixth curves bifurcating from the fourth hyperbola are members of $\mathcal{C}_6$ or $\mathcal{C}_1$. All other primary branches are similar to these examples. We can then follow these bifurcation curves and as they bifurcate to daughter branches with less symmetry, we can follow the path on the condensed symmetry digraph in an expected way. We know on what generation of bifurcation branches to find symmetry types in the digraph by how many edges away from the trivial symmetry type they are found. Figure~\ref{Strings}(a) shows a portion of the bifurcation diagram connecting the trivial branch through four bifurcations to a branch whose solution pairs are invariant under no symmetry transformations. The portion of (a) enclosed in the rectangle is enlarged in (b) for clarity, along with black diamonds marking the solutions whose contours are given in Figure~\ref{String} from right to left. \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}} \includegraphics[width=0.45\textwidth]{fig11a} % Stringsall.png &\includegraphics[width=0.45\textwidth]{fig11b} % StringsClose.png \\ (a) & (b) \end{tabular} \end{center} \caption{(a) is a part of the bifurcation diagram connecting the trivial branch to a branch whose solutions have no symmetry. The portion enclosed in the rectangle is enlarged in (b) for clarity. The black diamonds mark the nontrivial solutions whose contours are given in Figure~\ref{String} from right to left} \label{Strings} \end{figure} \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}} \includegraphics[width=0.14\textwidth]{fig12a} % String0.png &\raisebox{10mm}[0mm][0mm]{$\overset{\longrightarrow}{D_4}$} &\includegraphics[width=0.14\textwidth]{fig12b} % String1.png &\raisebox{10mm}[0mm][0mm]{$\overset{\longrightarrow}{Z_2}$} &\includegraphics[width=0.14\textwidth]{fig12c} % String2.png &\raisebox{10mm}[0mm][0mm]{$\overset{\longrightarrow}{Z_2}$} &\includegraphics[width=0.14\textwidth]{fig12d} % String3.png &\raisebox{10mm}[0mm][0mm]{$\overset{\longrightarrow}{Z_2}$} &\includegraphics[width=0.14\textwidth]{fig12e} % String4.png \\ $\langle \rho,\tau,\sigma,-1\rangle$ &&$\langle -\rho^2,\tau,\sigma\rangle$ &&$\langle -\rho^2,\tau\rangle$ &&$\langle -\rho^2\rangle$ &&$\langle 1\rangle$ \\ $\in\mathcal{C}_{0}$ &&$\in\mathcal{C}_{3}$ &&$\in\mathcal{C}_{10}$ &&$\in\mathcal{C}_{15}$ &&$\in\mathcal{C}_{20}$ \end{tabular} \end{center} \caption{The contour plots of solution pairs that correspond to the black diamonds in Figure~\ref{Strings}(b), from right to left, given along with their symmetry types and condensation classes.} \label{String} \end{figure} \begin{figure}[ht] \begin{center} \begin{tabular}{@{}c@{}c@{}c@{}c@{}c@{}} \includegraphics[width=0.37\textwidth]{fig13a} % JValCCN.png &\includegraphics[width=0.14\textwidth]{fig13b} % C3_S78.png &\includegraphics[width=0.14\textwidth]{fig13c} % C3_S60.png &\includegraphics[width=0.14\textwidth]{fig13d} % C3_S66.png &\includegraphics[width=0.14\textwidth]{fig13e} % C3_S77.png \\ (a) & (b) & (c) & (d) & (e) \end{tabular} \end{center} \caption{A plot of the value of $J_{\lambda,\mu}(\bm{U})$ for the four primary branches whose solutions change sign exactly once is given in (a). Figures (b) and (c) are the two solution pairs from (a) for parameter values of $\lambda=20$ from bottom to top. Figures (d) and (e) are the two remaining solution pairs for $\lambda=-60$ from bottom to top. The solution pairs shown in (b) and (d) appear to be the lowest energy among sign changing solutions, hence their existence might be proven borrowing techniques from~\cite{CCN}} \label{CCN} \end{figure} \section{Future direction} \label{Future} A proof of the existence of a minimal energy sign-changing exactly-once solution for this system is of interest. Due to the infinite number of negative eigenvalues when dealing with the full basis of $H\times H$, the concept of Morse index is not defined, meaning an alternate to the mountain pass theorem must be used. To motivate future effort in this direction, Figure~\ref{CCN} shows the value of $J_{\lambda,\mu}$ for the different bifurcation branches of the possible candidates for such a solution pair and the corresponding contour plots of the solutions on these branches. We find that the solution type of the lowest energy at $\lambda=20$ seems related to the one proven to exist in~\cite{CCN} for the single PDE case. We propose representing the Galerkin basis by a rotation of the chosen basis vectors that lends itself to the types of solutions present in the positive and negative parameter spaces, respectively. If we define $X=\textrm{span}\{(\psi_i,\psi_i)\}$ and $Y=\textrm{span}\{(\psi_i,-\psi_i)\}$. Then $H\times H=X\oplus Y$ and the Hessian of $J_{0,0}(\bm{0})$ will be negative definite on $X$ and positive definite on $Y$. In this way, we might be able to treat the two sets separately, each with a different sense of Morse index. \subsection*{Conclusions} We have presented a numerical method for studying Hamiltonian systems of BVP. These are most likely to arise as steady states of coupled systems of reaction diffusion equations with nonlinear reaction terms. In addition, we gave an analysis of the symmetry of the space of pairs of functions defined on the unit square, along with a representative example to illustrate the main results. Finally, we have proposed a future direction for this work to take a more theoretical approach in proving the existence of a minimal energy sign-changing exactly-once solution to the system where Morse index is not well defined due to infinitely many negative eigenvalues. \begin{thebibliography}{99} \bibitem{CCN} Alfonso Castro, Jorge Cossio, John~M. Neuberger; \newblock A sign-changing solution for a superlinear {D}irichlet problem. \newblock {\em Rocky Mountain J. Math.}, 27(4):1041--1053,1997. \bibitem{Defig} Djairo~G. de~Figueiredo; \newblock {\em Semilinear elliptic systems: existence, multiplicity, symmetry of solutions}. \newblock Handb. Differ. Equ. Elsevier/North-Holland, Amsterdam, 2008. \bibitem{GAP} The GAP~Group; \newblock {\em {GAP -- Groups, Algorithms, and Programming, Version 4.4.9}}, 2006. \bibitem{NSS3} John~M. Neuberger, N{\'a}ndor Sieben, James~W. Swift; \newblock Automated bifurcation analysis for nonlinear elliptic partial difference equations on graphs. \newblock {\em Internat. J. Bifur. Chaos Appl. Sci. Engrg.}, 19(8):2531--2556, 2009. \bibitem{NS} John~M. Neuberger, James~W. Swift. \newblock Newton's method and {M}orse index for semilinear elliptic {PDE}s. \newblock {\em Internat. J. Bifur. Chaos Appl. Sci. Engrg.}, 11(3):801--820, 2001. \end{thebibliography} \end{document}