The \((3+1)\)-dimensional evolution of an inhomogeneous axion field configuration around the QCD epoch is studied numerically, including important non-linear effects due to the attractive self-interaction. It is found that axion perturbations on scales corresponding to causally disconnected regions at \(T \sim 1 \, {\rm GeV}\) can lead to very dense pseudo-soliton configurations we call axitons. These configurations evolve to axion miniclusters with present density \(\rho_a \ga 10^{-8}\,{\rm g \, cm^{-3}}\). This is high enough for the collisional \(2a \rightarrow 2a\) process to lead to Bose--Einstein relaxation in the gravitationally bound clumps of axions, forming Bose stars.