## 5 Evolution  
While the solutions found with the GA serve as good approximations to the exact solutions, they are not exact. Let us denote the exact solution of the eigenvalue problem as \(\psi_{\text{exact}}\) and the numerical solution as \(\psi_{\text{exact}}+\delta\psi\), where \(\delta\psi\) represents the error associated with the numerical truncation error of using a discrete domain for its construction. Stability can be tested by monitoring the behavior of this error over time when the exact solution is used as the initial condition plus the perturbation. Specifically, we analyze the dynamics triggered by such perturbation. The stationary solution is deemed stable if the perturbation remains bounded during its evolution, and it is considered unstable if the perturbation grows over time. This error analysis is commonly used to test the convergence of numerical solutions of stable systems [6] and to test how the errors converge to zero while numerical resolution is increased. Thus, the evolution of the numerical solution of the initial value problem has an error that we show does not diverge (see, e.g., [6; 13]). But we would like to monitor the error when using the empirical density profile as initial condition and see how it behaves. For this, we programmed a code that evolved these initial conditions by solving the time-dependent system (1) and (2) using spherical symmetry. The solution took place on the same numerical domain \(D^{d}\) used to solve the eigenvalue problem. The code used the method of lines to solve GP Equation (1) with a fourth-order Runge-Kutta integrator and second-order accurate stencils for spatial derivatives. At the origin, the order parameter was extrapolated with a second-order accurate approximation. Simultaneously, at each intermediate step of the Runge-Kutta method, we solved Poisson Equation (2) outwards from the origin until \(r_{max}\). The diagnostics we used to monitor the growth of perturbations was the central density at the origin, and the results are shown in Figure 4 for the evolution of the solution to the eigenvalue problem and the empirical profile (24) as well as that of Formula (29) obtained from simulations in [20], for various combinations of \(\psi_{c}\) and \(g\). In the left column, we show the time series of the central density, while the right column shows its Fourier transform that illustrates the triggered oscillation modes. The first case \((\psi_{c},g)=(2,-0.5)\) is an unstable solution that collapses due to the attractive nature of self-interaction (negative \(g\)). In this case, the density when initial conditions are the numerical solution blows up at a finite time, whereas when evolving the empirical profile, the density also acquires an out-of-bound central density, also indicating  
Figure 3: Solutions at the vertices of the parameter space. The top left plot corresponds to \((\psi_{c},g)=(2,-0.5)\), the top right plot to \((\psi_{c},g)=(2,0.5)\), the bottom left plot to \((\psi_{c},g)=(1,-0.5)\), and the bottom right plot to \((\psi_{c},g)=(1,0.5)\). Points represent the solution to the eigenvalue problem, while continuous lines represent our empirical Formula (24) and dashed lines empirical Formula (29) obtained from simulations in [20]. instability. In that case the Fourier spectrum says little about the oscillation modes and is not shown. There are also three stable cases with \((\psi_{c},g)=(1,-0.5)\), \((2,\pm 0.5)\), whose central density oscillates with different amplitudes and frequency modes. As expected, the density of the solution of the eigenvalue problem is closer to the exact solution than empirical Formula (24). An implication is that in the first case, the amplitude of the oscillations triggered by the truncation error is smaller than in the second case, where the difference between the numerical solution of the eigenproblem and the empirical formula add an extra perturbation. The amplitudes differ approximately by an order of magnitude.