Perhaps a good way to understand this is through the example of the $H_2$ molecule. Before delving into the molecule, however, let us examine the Hydrogen atom more closely. The $1s$ state of the Hydrogen atom is filled by one electron, but it can take two at most because of spin. Now, the other Hydrogen atom also has one atom in the $1s$ state and can again accommodate at most two.
If we were to bring the two atoms closer together, and then ultimately form a molecule, do you think that the number of states changes? Of course not! The molecule again must be able to accommodate a total of four electrons (two from each atom). (We are ignoring higher energy states for now).
When we think about the Hydrogen molecule, there is a very strong constraint on what the wavefunction must look like. This is because the Hydrogen molecule possesses parity or reflection symmetry. This means that when I switch the atoms, I must either get a $+$ sign or a $-$ sign on the wavefunction. (See the figure below.)
(source: tu-bs.de)
To be more explicit, let's represent the $1s$ states on the Hydrogen atoms 1 and 2 when they are far apart by $|1s_{1}\rangle$ and $|1s_{2}\rangle$. When we bring the atoms close together, the reflection symmetry constraint imposes a severe constraint. Let's see how we can combine the wavefunctions of the individual atoms to satisfy this constraint. It turns out that the following combinations work quite well:
\begin{equation}
|\psi_B\rangle = |1s_1\rangle + |1s_2\rangle
\end{equation}
\begin{equation}
|\psi_A\rangle = |1s_1\rangle - |1s_2\rangle
\end{equation}
where $B$ and $A$ denote bonding and anti-bonding respectively. Now look what happens when I switch the places of the two molecules by applying the reflection operation, $R$:
\begin{equation}
R|\psi_B\rangle = |1s_2\rangle + |1s_1\rangle = +|\psi_B\rangle
\end{equation}
\begin{equation}
R|\psi_A\rangle = |1s_2\rangle - |1s_1\rangle = -|\psi_A\rangle
\end{equation}
This set of wavefunctions therefore satisfies the constraint that we set by imposing reflection symmetry. Now it turns out that because the reflection symmetry, the wavefunctions must look something like this:
Again, each orbital can accommodate two electrons for a total of four, exactly the number of states we started out with. As a rule of thumb, the energy of the orbital with the greater curvature is usually the more energetic one, which is the anti-bonding one in this case. Therefore, the bonding orbital will be filled by the two electrons.
I think that the way you are trying to visualize the problem is flawed. You should not think of the standing wave as a sine wave extending through space infinitely. It is really cut off by the potential of the hydrogen atom and looks like the red and green curves above for the individual atoms. When you bring them closer together and they mix you get the purple and blue curves. Hope this helps!
You're right on a lot of counts. The wavefunction of the system is indeed a function of the form
$$
\Psi=\Psi(\mathbf r_1,\mathbf r_2),
$$
and there's no separating the two, because of the cross term in the Schrödinger equation. This means that it is fundamentally impossible to ask for things like "the probability amplitude for electron 1", because that depends on the position of electron 2. So at least a priori you're in a huge pickle.
The way we solve this is, to a large extent, to try to pretend that this isn't an issue - and somewhat surprisingly, it tends to work! For example, it would be really nice if the electronic dynamics were just completely decoupled from each other:
$$
\Psi(\mathbf r_1,\mathbf r_2)=\psi_1(\mathbf r_1)\psi_2(\mathbf r_2),
$$
so you could have legitimate (independent) probability amplitudes for the position of each of the electrons, and so on. In practice this is not quite possible because the electron indistinguishability requires you to use an antisymmetric wavefunction:
$$
\Psi(\mathbf r_1,\mathbf r_2)=\frac{\psi_1(\mathbf r_1)\psi_2(\mathbf r_2)-\psi_2(\mathbf r_1)\psi_1(\mathbf r_2)}{\sqrt{2}}.
\tag1
$$
Suppose that the eigenfunction was actually of this form. What could you do to obtain this eigenstate? As a fist go, you can solve the independent hydrogenic problems and pretend that you're done, but you're missing the electron-electron repulsion. You could solve the hydrogenic problem for electron 1 and then put in its charge density for electron 2 and solve its single electron Schrödinger equation, but then you'd need to go back to electron 1 with your $\psi_2$. You can then try and repeat this procedure for a long time and see if you get something sensible.
Alternatively, you could try reasonable guesses for $\psi_1$ and $\psi_2$ with some variable parameters, and then try and find the minimum of $⟨\Psi|H|\Psi⟩$ over those parameters, in the hope that this minimum will get you relatively close to the ground state.
These, and similar, are the core of the Hartree-Fock methods. They make the fundamental assumption that the electronic wavefunction is as separable as it can be - a single Slater determinant, as in equation $(1)$ - and try to make that work as well as possible. Somewhat surprisingly, perhaps, this can be really quite close for many intents and purposes. (In other situations, of course, it can fail catastrophically!)
In reality, of course, there's a lot more to take into account. For one, Hartree-Fock approximations generally don't account for 'electron correlation' which is a fuzzy term but essentially refers to terms of the form $⟨\psi_1\otimes\psi_2| r_{12}^{-1} |\psi_2\otimes\psi_1⟩$. More importantly, there is no guarantee that the system will be in a single configuration (i.e. a single Slater determinant), and in general your eigenstate could be a nontrivial superposition of many different configurations. This is a particular worry in molecules, but it's also required for a quantitatively correct description of atoms.
If you want to go down that route, it's called quantum chemistry, and it is a huge field. In general, the name of the game is to find a basis of one-electron orbitals which will be nice to work with, and then get to work intensively by numerically diagonalizing the many-electron hamiltonian in that basis, with a multitude of methods to deal with multi-configuration effects. As the size of the basis increases (and potentially as you increase the 'amount of correlation' you include), the eigenstates / eigenenergies should converge to the true values.
Having said that, configurations like $(1)$ are still very useful ingredients of quantitative descriptions, and in general each eigenstate will be dominated by a single configuration. This is the sort of thing we mean when we say things like
the lithium ground state has two electrons in the 1s shell and one in the 2s shell
which more practically says that there exist wavefunctions $\psi_{1s}$ and $\psi_{2s}$ such that (once you account for spin) the corresponding Slater determinant is a good approximation to the true eigenstate. This is what makes the shells and the hydrogenic-style orbitals useful in a many-electron setting.
However, a word to the wise: orbitals are completely fictional concepts. That is, they are unphysical and they are completely inaccessible to any possible measurement. (Instead, it is only the full $N$-electron wavefunction that is available to experiment.)
To see this, consider the state $(1)$ and transform it by substituting the wavefunctions $\psi_j$ by $\psi_1\pm\psi_2$:
\begin{align}
\Psi'(\mathbf r_1,\mathbf r_2)
&=\frac{\psi_1'(\mathbf r_1)\psi_2'(\mathbf r_2)-\psi_2'(\mathbf r_1)\psi_1'(\mathbf r_2)}{\sqrt{2}}
\\&=\frac{
(\psi_1(\mathbf r_1)-\psi_2(\mathbf r_1))(\psi_1(\mathbf r_2)+\psi_2(\mathbf r_2))
-(\psi_1(\mathbf r_1)+\psi_2(\mathbf r_1))(\psi_1(\mathbf r_2)-\psi_2(\mathbf r_2))
}{2\sqrt{2}}
\\&=\frac{\psi_1(\mathbf r_1)\psi_2(\mathbf r_2)-\psi_2(\mathbf r_1)\psi_1(\mathbf r_2)}{\sqrt{2}}
\\&=\Psi(\mathbf r_1,\mathbf r_2).
\end{align}
That is, the Slater determinant that comes from linear combinations of the $\psi_j$ is indistinguishable from the one you get from the $\psi_j$ themselves. This extends to any basis change on that subspace with unit determinant; for more details see this thread. The implication is that labels like s, p, d, f, and so on are useful to describe the basis functions that we use to build the dominating configuration in a state, but they cannot be reliably inferred from the many-electron wavefunction itself. (This is as opposed to term symbols, which describe the global angular momentum characteristics of the eigenstate, and which can indeed be obtained from the many-electron eigenfunction.)
Best Answer
It means that the theoretical explanation that you are being given in chemistry is incomplete. I wouldn't call it false because it is a very useful approximation that explains pretty much everything that we need in entry level chemistry "well enough" (even quantitatively), but it's just not how quantum mechanics of many-particle systems actually works. Let's try to clarify what is really happening.
The orbitals that you are used to from the drawings in physics and chemistry books are the single electron wave functions of hydrogen. In classical physics these orbitals would indeed "overlap" for a multi-electron wave function. Quantum mechanics doesn't work that way, though. A multi-particle wave function doesn't live in the same Hilbert space as a single particle wave function. Instead it lives in the product of several such Hilbert spaces, i.e. it is a higher dimensional distribution than the single electron function. For a single electron atom like hydrogen we can write the wave function as ψ(r). For a two electron atom, however, it becomes ψ(r1,r2), i.e. it now depends on six coordinates rather than three. For n-electrons we have to use ψ(r1,r2,⋯,ri), so we need a 3n-dimensional function to represent the system correctly.
That we can draw orbitals for the hydrogen problem is a mere coincidence because the spatial dimensionality of the solution happens to have the same three dimensions as space itself. For two electrons, however, the actual solution of the corresponding Schroedinger equation would live in six dimensions and so on. This can not be visualized easily. In the "orbital approximation" we are replacing this higher dimensional solution with a much more simple product:
ψ(r1,r2,⋯,ri)≈φ1(r1)φ2(r2)⋯φi(ri)
and we further assume that the electrons (described by the φi(ri) functions) behave independently like in a single electron atom. The reason why we are doing this is because the Schroedinger equation for the full 3n dimensional many-body problem can't be solved in closed form. Even the approximation problem is hard, but it can be analyzed sufficiently well to extract important information about atoms other than hydrogen and even molecules.