## 3 Numerical Method and Solutions
### Discretization Scheme  
We outline the process to numerically find equilibrium solutions of System (2). We take advantage of working with a bounded, periodic grid and use the Fast Fourier Transform to move into frequency space, compute derivatives and convolutions, and use the Inverse Fast Fourier Transform to move back to physical space. This process gives us a discrete representation of the right hand side of our system, which allows us to discretize in time and numerically solve the system. Consider the Fourier Transform in two dimensions,  
\[\hat{a}(k_{1},k_{1})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}u(x_{1},x_{ 2})e^{-2\pi i(k_{1}x_{1}+k_{2}x_{2})}dxdy,\ \ k_{1},k_{2}\in\mathbb{R}.\]  
Because we are solving a nonlocal system of PDEs, we must address computing convolutions and derivatives efficiently. Therefore, there we employ two convenient properties of the Fourier Transform in our computations:  
\[\widehat{\partial_{x_{i}}f}(x_{1},x_{2}) =ik_{i}\widehat{f}(k_{1},k_{2}) \tag{4}\] \[\widehat{K\ast f}(k_{1},k_{2}) =\widehat{K}(k_{1},k_{2})\widehat{f}(k_{1},k_{2}). \tag{5}\] \begin{table}
\begin{tabular}{c c} \hline
**Term in Equation** (2) & **Physical Description** \\ \hline \(\partial_{t}u_{i}(x,t)\) & Rate of change in population density \(i\) with time \\ \hline \(\eta\Delta u_{i}^{2}(x,t)\) & Diffusion modeling an intra-group overcrowding effect \\ \hline \(u_{i}\
abla(K_{1}\ast u_{i})\) & Intra-group non-local aggregation \\ \hline \(-u_{i}\
abla(K_{2}\ast u_{i})\) & Inter-group non-local segregation \\ \hline \(u_{i}\
abla U(x)\) & Environmental advection \\ \hline \end{tabular}
\end{table} 
Table 1: Physical description of the terms in the model given by Equation (2). As we are working on a bounded, periodic, \(M\times M\) grid, we can discretize in \(x\)\(y\),and \(t\) such that \(x_{m}=m\Delta x\), \(y_{j}=j\Delta y\), where \(m,j\in\{0,\ldots...,M-1\}\), and \(t_{n}=n\Delta t\), \(n\in\mathbb{N}\). Thus, \(u(x_{m},y_{j},t_{n})=u_{mj}^{n}\). Below, we drop the dependency on \(n\) for notational convenience. For a fixed time \(t_{n}\), we compute derivatives and convolutions in frequency space using the Discrete Fourier Transform,  
\[\hat{u}_{k,v}=\frac{1}{M^{2}}\sum_{m=0}^{M-1}\sum_{j=0}^{M-1}u_{mj}e^{-\frac{2 ni}{M}(km+v)j},\ \ k,v\in\{0,...,M-1\}, \tag{6}\]  
which allows us to pass into frequency space. The Inverse Discrete Fourier Transform,  
\[u_{mj}=\sum_{k=0}^{M-1}\sum_{v=0}^{M-1}\hat{u}_{k,v}e^{\frac{2 ni}{M}(km+v)j},\ \ m,j\in\{0,1,...,M\}, \tag{7}\]  
allows us to return to physical space. We utilize Equation (6) to move into the Fourier domain, we use properties (4) and (5) to compute derivatives and convolutions in Fourier space, and use (7) to move back to physical space. In fact, we make use of the Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT) to move in and out of frequency space efficiently, reducing number of operations from \(M^{2}\) to \(M\log M\). We first compute the right-hand side of System (2) on the bounded, periodic grid for a fixed \(t=n\Delta t\). Let \(u_{mj}^{i,m}\) be the value of \(u_{i}\) at time \(n\Delta t\) for group \(i\in\{1,...,N\}\) at grid point \((m,j)\) for \(m,j\in\{0,1,...,M-1\}\), and let \(\hat{u}_{kr}^{i,n}\) represent the DFT for \(k,v\in\{0,1,...,M-1\}\). We begin by using the FFT to move into Fourier space to determine \(\hat{u}_{kv}^{i,n}\), \(\hat{K}\), and \(\hat{U}\). We compute the nonlocal terms in System (2) by multiplying \(\hat{u}_{kj}^{i,n}\) and \(\hat{K}\) for each group \(i\), utilizing the property in Equation (5). We remain in Fourier space to compute the gradient of the nonlocal terms and \(U\), using the property in Equation (4),  
\[\left(\widehat{\
abla K*u}\right)_{kv}^{n}=(ik\hat{K}_{kv}\hat{u}_{kv}^{n},iv \hat{K}_{kv}\hat{u}_{kv}^{n}),\]  
and  
\[\left(\widehat{\
abla U}\right)_{kv}^{n}=(ik\hat{U}_{kv}^{n},iv\hat{U}_{kv}^{ n}).\]  
We use the IFFT to return the physical domain and multiply these terms by \(u_{mj}^{i,n}\).