I will provide an answer in multiple parts, and most of the information I describe here comes from source [2].
For a $\ce{H2}$ molecule, the Hamiltonian must consider the interactions between both electrons with both nuclei. So the condensed version would be as described in source [1] and applying the Born-Oppenheimer approximation:
$$\hat{H} = -\frac{1}{2}(\nabla^2_1 + \nabla^2_2) - \frac{1}{r_{1A}} - \frac{1}{r_{1B}} - \frac{1}{r_{2A}} - \frac{1}{r_{2B}} + \frac{1}{r_{12}} + \frac{1}{R}$$
Here, $r_{iA}$ represents how far electron $i$ is from atom A and atom B; $i=1,2$. $r_{12}$ accounts for the inter-electronic interaction and $R$ represents the distance between nuclei. While $R$ can be fixed, the energies are typically plotted against the internuclear distance like a function $E(R)$, hence where you get the Morse potential from.
History and Methods
We will focus on the spatial part for now. In 1927, Heitler and London proposed a method to analyze $\ce{H2}$ in particular called the valence-bond method. The point here is that we'll be constructing a trial function to approximate the actual wave function. There are several different approaches, one of which involves the linear combination of atomic orbitals with molecular-orbital method, but I'll get into that later. First, let's consider a trial wavefunction between two hydrogen atoms that are infinitely apart:
$$\psi_1 = 1s_A(1)1s_B(2)$$
Where $1s_A(1)$ represents electron 1 on the 1s orbital of atom A and $1s_B(2)$ represents electron 2 on the 1s orbital of atom B. One thing we have to consider here is that when the atoms are put closer together, the electrons are indistinguishable. That means we must also consider
$$\psi_2 = 1s_A(2)1s_B(1)$$
For example, electron 2 may be on $1s_A$ and electron 1 on $1s_B$ from the viewpoint of atomic orbitals. This means that the overall wave function can be represented as a linear combination of the two
$$\Psi_{VB} = c_1\psi_1 + c_2\psi_2$$
However, we should also include an ionic term. For example, if two electrons are on one side of the molecule, you'll have a positively and negatively charged molecules on either sides.
$$\Psi_{ionic} = 1s_A(1)1s_A(2) + 1s_B(1)1s_B(2)$$
We will come back to this, but we need to first talk about $\ce{H2+}$ and molecular orbitals.
$\ce{H2+}$ Molecular Orbitals
In molecular-orbital method, we essentially derive the molecular wave functions as a product of molecular orbitals (which are single-electron molecular wave function, see source [4]).
So, to understand how we came to the trial wave function for $\ce{H2}$, we need to first look at $\ce{H2+}$. The Hamiltonian is
$$\hat{H} = -\frac{1}{2}\nabla^2 - \frac{1}{r_A} - \frac{1}{r_A} + \frac{1}{R}$$
Where we have two nuclei A and B and one electron. While it is possible to solve this system exactly, we will solve it approximately to construct the molecular wave function for $\ce{H2}$ eventually. So take:
$$\psi_{\ce{H+}} = c_11s_A + c_21s_B$$
This is a linear combination of atomic orbitals. Again, the point is to create a trial wave function that can approximate the actual wave function of $\ce{H2+}$ well. Taking a linear combination of atomic orbitals has shown success, as we will later see in the treatment of $\ce{H2}$, and that is why it is widely used. The reason why we take the linear combination from both atoms is because they both contribute to the overall wave function.
For this trial wave function, the associated secular determinant would be
$$\begin{vmatrix}
H_{AA}-E\cdot1 & H_{AB} - RS\\
H_{AB}-ES & H_{BB}-E\cdot1\\
\end{vmatrix}$$
When we solve this, though it is somewhat complex and I don't want to type it all out, we will get two solutions, corresponding to the molecular orbitals, bonding and antibonding, formed with an $\ce{H2+}$ molecule. They are
$$\psi_+ = \frac{1}{\sqrt{2(1+S)}}(1s_A + 1s_B)$$
$$\psi_- = \frac{1}{\sqrt{2(1+S)}}(1s_A - 1s_B)$$
where $S = \int dr1s_A1s_B$. Source [2] provides the pathway to the solutions. Here, we've developed a set of molecular orbitals that are analogous to $\ce{H2}$, and we can use them to tackle it.
$\ce{H2}$ Analysis
What the molecular-orbital method wants us to consider is constructing a molecular wave function by combining molecular orbitals. With $\ce{H2+}$, we found the molecular orbitals using a linear combination of atomic orbitals as a trial function. Since $\psi_+$ is the molecular orbital with the ground state energy, we can basically describe $\ce{H2}$ by filling that bonding molecular orbital with another electron with opposite spin. This can be represented by a Slater determinant
$$\Psi = \begin{vmatrix}
\psi_+\alpha(1) & \psi_+\beta(1)\\
\psi_+\alpha(2) & \psi_+\beta(2)\\
\end{vmatrix}$$
$$= \psi_+\cdot\psi_+[\frac{1}{\sqrt{2}}[\alpha(1)\beta(2)-\alpha(2)\beta(1)]]$$
Considering only the spatial part, we get
$$\Psi_{MO} = \frac{1}{2(1+S)}[1s_A(1) + 1s_B(1)][1s_A(2) + 1s_B(2)]$$
$$ = 1s_A(1)1s_B(2) + 1s_A(2)1s_B(1) + 1s_A(1)1s_A(2) + 1s_B(1)1s_B(2)$$
Interestingly, if you do the math, this is almost just like combining the valance-bond method with the ionic terms but without the coefficients. This method is called the LCAO-MO method. This essentially explains how we got to molecular orbital theory, why we express molecular orbitals as a linear combination (as in the section about $\ce{H2+}$), and how we use that to analyze $\ce{H2}$.
Although this is the basis of molecular orbital theory, it is the simple version. We should technically extend molecular-orbital theory as the in our earlier equation with $\Psi_{VB}$ and $\Psi_{ionic}$, each term isn't contributing equally, so we should add coefficients to each term in our trial wave function and solve variationally. This is done by considering the excited antibonding state $\psi_-$, and, amazingly, when we do, we get precisely the form $\Psi_{CI} = c_1\Psi_{VB} + c_2\Psi_{ionic}$. The subscript in $\Psi_{CI}$ represents the extension of simple molecular-orbital theory by including excited-state configurations, known as Configuration Interaction.
In all, we could have used $\Psi_{VB}$, $\Psi_{MO}$ (simple version), and $\Psi_{CI}$ (extended version), or any other more complicated trial wave function to analyze $\ce{H2}$. It's just that the molecular-orbital method with configuration interactions considered is widely used. Source [2] uses this extended form to analyze $\ce{H2}$, but it's rather complicated.
Note that we are not taking into account the $\ce{2s}$ orbital in this analysis as are considering the minimal basis set which leads to two molecular orbitals. Even in molecular orbital diagrams, you will see that the $\ce{H2}$ molecular orbitals arise from a linear combination of the 1s orbitals. We could choose a different trial function to approximate the wavefunction with for example: $\Psi = c_11s_A + c_21s_B + c_32s_A + c_42s_B$. What we find is that when variationally computed, the ground-state molecular orbital becomes $\Psi = 0.7071(1s_A + 1s_B) + 0.00145(2s_A + 2s_B)$. The corresponding energy basically has no improvement compared to solely analyzing the linear combinations of $\ce{1s}$ orbitals with the valance-bond method (see source [2]). So, we are consider solely $\ce{1s}$ here.
Your Analysis
This derivation is in ways similar and different than yours. You basically have the form
$$\Psi_T = [c_11s_A(1) + c_21s_B(2)][c_31s_B(1) + c_41s_B(2)]$$
Sure, you could go about solving this by partitioning the Hamiltonian and ignoring the electron-electron repusion term, but all you really get is the $\ce{H2+}$ wavefunction for each electron. Also remember that the energies depend on the internuclear distance, and that is considered in the analysis of $\ce{H2+}$. So if we also don't include that in the Hamiltonian, then we aren't really analyzing a molecule anymore. Plus, the hard part is not constructing the secular determinant, the hard part is try to mathematically analyze the integrals inside the secular determinant.
As for your comment on spin, for this investigation, we say that the Hamiltonian does not depend on the spin, so the energy is determined by the spatial elements of the wave function.
Sources
[1] http://websites.umich.edu/~chem461/QMChap10.pdf
[2] Quantum Chemistry - Donald A. McQuarrie
[3] https://aip.scitation.org/doi/pdf/10.1063/1.1701472?casa_token=If9tKkNIyIQAAAAA:xYrx4IoQExf2NFTRSNA-ElywJRrhTmbsRrYZIlSPx_vZS9e9DI5EHcMKyU2jJ0AqrIU2d821bbE
[4]http://lampx.tugraz.at/~hadley/ss1/molecules/mo.php#:~:text=In%20the%20simplest%20approximation%2C%20the,r%E2%88%92%E2%83%97ra%7C.
Other Sources
[5] https://www.sciencedirect.com/topics/mathematics/london-theory#:~:text=The%20Heitler%E2%80%93London%20(Ref.,as%20shown%20in%20Figure%2013.2.
[6] https://dasher.wustl.edu/chem478/reading/frontier-orbitals-anh.pdf
[7] http://websites.umich.edu/~chem461/QMChap10.pdf
[8]https://chem.libretexts.org/Courses/University_of_California_Davis/UCD_Chem_107B%3A_Physical_Chemistry_for_Life_Scientists/Chapters/5%3A_The_Chemical_Bond/5.2%3A_Valence_Bond_Theory#:~:text=along%20the%20axis.-,The%20Valence%20Bond%20Wavefunction,Equation%205.2.
[9] https://www.sciencedirect.com/topics/chemistry/molecular-orbital-method
[10] https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1951.0240
