## 4 Three-Dimensional Numerical Results  
In this section, we consider the full 3D model for direct numerical simulations of Equation (20). First, we find the ground-state solution by solving Equation (20) by means of the imaginary-time (\(t\rightarrow-\mathrm{i}t\)) propagation in the presence of the external trap and two-body contact plus dipolar interactions in the absence of three-body interactions. Then, the evolution of the input provided by the ground state is simulated in the real-time. In the course of the real-time evolution, the external trap is adiabatically removed and the three-body interaction is ramped up which brings the system to the same form for which the VA was elaborated above. We now consider the following 3D GP equation  
\[ \begin{split} i\frac{\partial\phi(\mathbf{r},t)}{\partial t}=& \bigg{[}-\frac{1}{2}\
abla^{2}+d(t)\frac{1}{2}(\gamma^{2}x^{2}+ \
u^{2}y^{2}+\lambda^{2}z^{2})+g(t)|\phi(\mathbf{r},t)_{i}|^{2}+\\ &\chi(t)|\phi(\mathbf{r},t)|^{4}+3 Na_{\mathrm{dd}}\int\frac{1-3 \cos^{2}\theta}{|\mathbf{R}|^{3}}|\phi(\mathbf{r}^{\prime},t)|^{2}d\mathbf{r }^{\prime}\bigg{]}\phi(\mathbf{r},t).\end{split}  \tag{20}\]  
and numerically solve by gradually decreasing \(d(t)\) and increasing \(g(t)\), \(\chi(t)\). The DD interaction was evaluated by means of the fast Fourier transform [8]. Typical space and time steps for the numerical grid were taken as 0.05 and 0.005, respectively. In the course of time evolution, we increase the nonlinearity coefficients from 0 at each time step as \(g(t)=f(t)g_{f}\{a_{1}-b_{1}\sin[\Omega\left(t-\tau\right)]\}\) and \(\chi(t)=f(t)\chi_{f}\{a_{2}-b_{2}\sin[\Omega\left(t-\tau\right)]\}\), with \(f(t)=t/\tau\) for \(0\leq t\leq\tau\) and \(f(t)=1\) for \(t>\tau\). At the same time, the trap is removed by changing \(d(t)\) from 1 to 0 by \(d(t)=1-f(t)\). Thus, the trap is eventually removed and the final values of the coefficients are attained at \(t=\tau\) and the periodically oscillating nonlinearities become \(g(t)=g_{f}[a_{1}-b_{1}\sin(\Omega t)]\) and \(\chi(t)=\chi_{f}[a_{2}-b_{2}\sin(\Omega t)]\) is introduced for \(t>\tau\), cf. Refs. [33, 34, 36]. After the complete removal of the external trap, the condensates will collapse for strong attraction due to the final nonlinearity. If the attraction due to the final value of nonlinearity after switching off the external trap is too weak, the system becomes weakly attractive and the condensates start to expand. In Figure 6, we present the results for the condensate of \({}^{52}\)Cr atoms which have a moderate dipole moment corresponding to \(a_{dd}=16a_{0}\)[1, 2]. The figure demonstrates how one can stabilize the condensate with \(\epsilon_{0}=-25a_{0}\) by varying the strengths of the repulsive three-body interactions (without the contribution of the oscillatory part of the interactions). In panel (1,1) the condensates are found to be stable up to \(2 ms\) for \(\chi_{0}=-0.013\epsilon_{0}\). Similarly, in panels (1,2), (2,1) and (2,2), the condensates are stable up to \(2.5 ms\), \(3 ms\) and \(5 ms\) for \(\chi_{0}=-0.012\epsilon_{0}\), \(-0.011\epsilon_{0}\) and \(-0.01085\epsilon_{0}\), respectively. Further evolution of the condensates may perhaps lead to their collapse or expansion.