In a diatomic molecule there are four quantum mechanical particles. This is a very complex system and unlike for the single electron Schroedinger equation (aka the hydrogen problem) there are no known exact solutions. This is similar to the n-body problem in classical mechanics, which is also not solvable in closed form for n>2. To work around this difficulty we make several approximations to the full equations of motion.
The first one is called the Born-Oppenheimer approximation, in which we treat the movement of the nuclei in an effective field created by the electrons. This is possible because the nuclei are thousands of times heavier than the electrons and their movements are slow. The rotational and vibrational energy levels created by these movements are much smaller than the molecular binding energy.
Next we make a independent electron approximation, in which we treat each electron as moving in an effective field created by the nuclei and the other electron. This reduces the number of actual quantum mechanical states greatly. In mathematical terms we assume that we can factorize the multi-particle wave function $\psi(r_1, r_2)$ into a product of the form $\psi_1(r_1)\psi_2(r_2)$. This allows us to separate the many-particle Schroedinger equation into independent equations for each particle. Because we have to treat electrons as indistinguishable particles we actually end up with symmetrized and anti-symmetrized wave functions.
Finally, since this is still too hard to solve, we use a perturbation expansion of atomic orbitals (for which we have explicit solutions) and assume that we can combine them into molecular orbitals and that's how we estimate the energy levels.
So why do so many approximations work? Aren't we allowing for significant errors with this calculation? Yes, we are, however, the Born-Oppenheimer approximation is well motivated by the time-scale and energy separation between rotational vibration and electronic energy states.
The separation of the electronic wave function is possible because we are only interested in the lowest and most stable orbitals, i.e. those which are formed from e.g. two s-states and we are willing to neglect the Coulomb interaction terms between the two electrons. The s-p, p-p, s-d, p-d etc. molecular orbital states have higher energies (at least usually) and they are short lived excited states and since the electrons are farther away from the nuclei on average, their interaction terms are more important than for the s-s case. If we were interested in molecular fluorescence or reaction-dynamics, we could neither omit the interaction term nor could we neglect the higher orbitals in our considerations.
The symmetrization and anti-symmetrization is actually forced on us by the rules of quantum mechanics. The final approximation, to use perturbation series of atomic wave functions to describe molecular ones is the least well supported and there are extension of the method which use better wave functions to achieve a higher level or precision in estimating the binding energy of molecules.