As an experimentalist, I might not be the person best suited to answer this, but I'll give it a try.
The wavefunction is going to be difficult to visualize because in general it is a complex function. If you want to 'see' sqrt(-1), I suggest you resort to drugs, lots and lots of drugs. But as for physical interpretation, Born tells us that the amplitude of the wavefunction squared will give you the real probability distribution. So the wavefunction represents a probability amplitude, and the random nature arises in the measurement process.
When making a BEC, you extract energy from the atoms in a few ways. Usually, you start with laser cooling (in a magneto-optical trap, or MOT) that captures and cools about 1e9 atoms to about 100 microkelvin. From there, you often need to employ other tricks to beat the so-called doppler limit set by the linewidth or your laser. Another limit is the "recoil limit" that says you will never cool an atom to have less energy than the photon it (randomly) emits. So you turn the lasers off, and trap these very cold atoms in other ways only to cool them further. Either magnetic fields, or far-off-resonance, focused lasers (optical traps) are used to capture the atoms, but now we are not scattering lots of light off them like we were in the MOT. In both cases you are altering the energy landscape experienced by the cloud, and this causes it to seek the lowest energy. Lower the trap walls, and the hottest (fastest) atoms will escape, leaving the remaining (slower) atoms to thermalize to a lower temperature. This is called evaporative cooling. In magnetic traps, it is often done with radio frequencies which couple your atoms to different Zeeman sublevels, changing their potential energy, and ejecting them from the scene. In far detuned optical traps, they just lower the laser power.
As the atoms get colder, something remarkable happens. Their wave nature starts to emerge. The de Broglie wavelength gets longer as you get colder, until eventually, particle wavefunctions start overlapping with their neighbors. At this point, the atomic statistics are no longer well described by the good-ole Maxwell-Boltzmann distro, and we have to resort to Bose-Einstein statistics. (This of course assumes you are working with bosons, or integer spin particles.) What does that mean? Basically it means that particles start to favor occupying the same state, or 'mode'.
To determine the wavefunction of a BEC, you could start with the Schroedinger equation, where now you simply replace the single particle wave function with a many particle wavefunction that is built as a tensor product of the individual particles. This turns out to be a poor approach because the particles interact. They are constantly bumping into each other, and this manifests as a non-linear (i.e. density dependent) interaction term in the Hamiltonian. So we can far better approximate the situation with the Gross-Pitaevskii equation. If you want to model a BEC, the GPE might be a good place to start.
One of the striking features of QM in my mind is that the description of a particle is inextricably linked to its environment. It's right there in the GPE via the potential energy term, V. The solutions of this equation will be your wavefunction. When you confine atoms in a harmonic trap, one that looks like a parabola, your wavefunction will reflect this. Indeed it is not a gaussian you see, but an inverted parabola (see Thomas-Fermi approx.). It will look different in different traps. The residual gaussian is composed of thermal atoms, still described by the M-B picture.
I should note that when atoms are in the ground state, they still cannot be said to have "no motion". The uncertainty principle tells us that they will always smear out over some tiny volume of phase space. This is called the zero point motion.
The phase is a tricky subject. You are right that the book seems to have described it poorly. BECs are often said to be phase coherent, which is to say across the whole cloud, the phase will be the same. If it were not so, you would necessarily see striations since some of the atoms would destructively interfere with their neighbors. This is basically an excitation, and so energetically unfavorable. I might have just made a theorist lose her lunch with that explanation, but this whole thing is just meant to be a low level intro.
Good luck, and keep asking questions.
For a BEC, you want atoms to be in the same quantum state, not necessarily at the same position.
For a BEC, the temperature is low enough so that the de Broglie wavelength $\lambda_{\mathrm{dB}} \propto 1/\sqrt{T}$ is larger than the interatomic spacing $\propto n^{-1/3}$, $n$ being the density. This means that the wave nature of the atoms is large enough for it to be felt by other atoms, in other words atoms "see" each other even without exactly sitting on top of each other. This is just to further justify the claim that you don't need atoms at the same position. Actually, if you had a perfect box potential of side $L$, and you reached BEC, then the atoms will macroscopically occupy the ground state $ |\Psi|^2 \propto \sin^2(x/L)$ which is very much extended. If you let $L\rightarrow \infty$, the atomic distribution becomes flat. So, again, very much atoms not at the same positions.
Ok, so now interactions and collapse.
First of all, BEC is a non-interacting effect. It is not driven by a competition of interaction terms, but solely by Bose-Einstein statistics. It is experimentally interesting that BEC seems to exist also in interacting systems, though there is no general theoretical proof. By BEC in an interacting system I mean macroscopic occupation of the ground state + Off-Diagonal Long-Range Order (ODLRO) — so not all superfluids are BECs. Let me also point out that you need interactions to reach a BEC as you need to reach thermal equilibrium.
The interaction strength among weakly interacting Bose-condensed bosons is quantified by a $g n$ term in the Hamiltonian, where $g$ is $4\pi\hbar^2 a/m$ (Gross-Pitaevski equation). You can make this interaction attractive with $a<0$ and repulsive with $a>0$, where $a$ is the scattering length and it is given by $a(B) = a_0 f(B)$, where $a_0$ is the background scattering length in the presence on no external magnetic field $B$ ($f$ is some function).
The pressure of a weakly interacting Bose-condensed gas is (at $T=0$):
$$ P = -\frac{\partial E}{\partial V} = \frac{1}{2}gn^2.$$
Because $n^2$ is always positive, the condition for stability (i.e. not to collapse) is $P>0$ and hence $g>0 \Rightarrow a>0$ i.e. a repulsive system. With a positive pressure, the gas expands until it hits a wall (e.g. the confining potential). But if $P<0$ then the system is intrinsically unstable and collapses.
Rb-87 is "easy" because its background scattering length is positive and therefore trivially allows for a stable BEC. K-39, on the other hand, has a negative background scattering length so its "BEC" would collapse (and eventually explode). But its scattering length can be made repulsive by the use of a Feshbach resonance (applying a field $B$ to change $a$) so that it can undergo BEC.
Best Answer
As for your specific question about uncertainties, let's just do Galilean relativity: The observer at rest w.r.t. to the BEC has the momentum operator $p_0$ and another observer has a momentum operator $p_v$ and the difference $c = p_v - p_0$ is a constant (velocity of the observer times mass of the system under consideration). We have \begin{align} \sigma_{p_v} & = \sqrt{\langle p_v^2\rangle_- \langle p_v\rangle^2} = \sqrt{\langle (p_0 + c)^2\rangle_- \langle p_0 + c\rangle^2} \\ & = \sqrt{\langle p_0^2\rangle + 2\langle p_0\rangle c + \langle c^2\rangle - \langle p_0\rangle^2 - 2\langle p_0\rangle c - \langle c\rangle^2} \\ & = \sqrt{\langle p_0^2\rangle - \langle p_0\rangle^2 } = \sigma_{p_0} \end{align} that is, adding a constant to a random variable does not change its standard deviation, so at least in the Galilean case the position and momentum uncertainties of any given state do not depend on the observer, and so if you think of temperature as a measure of the momentum uncertainty of a substance, all Galilean observers agree on the temperature of a system.
In general the interaction of thermodynamics with special relativity is a bit controversial, see this question about the relativity of temperature in general and whether different observers might agree or disagree about the temperature of a system.