The differences between integrals evaluations of contracted an uncontracted basis functions
A lot of stuff is going on in this question. So for clarity let's make this answer very explicit. We can try to consider the two cases of using an uncontracted basis and a contracted basis to clarify the differences. As a case let's look at STO-$2$G for $\mathrm{H_2}$. In the STO-$2$G basis set $\mathrm{H_2}$ would normally be described by two basis functions each consisting of two primitives. 

Uncontracted case

If we use the STO-$2$G basis set in the case of using uncontracted Gaussian basis functions, each hydrogen will be described by two basis functions. Our overlap matrix will therefore be of the form:

$$ S= 
\begin{bmatrix}
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_1^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_2^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_3^\mathrm{H_2}\right.\right&gt;  &amp;
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_4^\mathrm{H_2}\right.\right&gt; \\
\left&lt;\phi_2^\mathrm{H_1}\left|\phi_1^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_2^\mathrm{H_1}\left|\phi_2^\mathrm{H_1}\right.\right&gt;  &amp;
\left&lt;\phi_2^\mathrm{H_1}\left|\phi_3^\mathrm{H_2}\right.\right&gt;  &amp; 
\left&lt;\phi_2^\mathrm{H_1}\left|\phi_4^\mathrm{H_2}\right.\right&gt;  \\
\left&lt;\phi_3^\mathrm{H_2}\left|\phi_1^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_3^\mathrm{H_2}\left|\phi_2^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_3^\mathrm{H_2}\left|\phi_3^\mathrm{H_2}\right.\right&gt;  &amp; 
\left&lt;\phi_3^\mathrm{H_2}\left|\phi_4^\mathrm{H_2}\right.\right&gt;  \\
\left&lt;\phi_4^\mathrm{H_2}\left|\phi_1^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_4^\mathrm{H_2}\left|\phi_2^\mathrm{H_1}\right.\right&gt;  &amp;
\left&lt;\phi_4^\mathrm{H_2}\left|\phi_3^\mathrm{H_2}\right.\right&gt;  &amp;
\left&lt;\phi_4^\mathrm{H_2}\left|\phi_4^\mathrm{H_2}\right.\right&gt; 
\end{bmatrix}
$$

In this notation we have that:

$$\left&lt;\phi_i^\mathrm{H_k}\left|\phi_j^\mathrm{H_l}\right.\right&gt; = 
N_iN_j\int_{-\infty}^\infty P_i\exp\left(-\alpha_ir_i^2\right)P_j\exp\left(-\alpha_jr_j^2\right) d\tau$$

With $P$ being a function of the angular momentum, in cartesian coordinates:

$$P=x^ly^mz^n$$

It can be seen that $P=1$ for s-functions. As we can see we have 4 basis functions in total (two on each hydrogen). You wrote something about, you would expect that in the case of uncontracted that we would have ""[number of electrons]x[number of basis functions] matrix"". As we can see this is not the case. Right now our matrix is in the atomic basis. When we solve the Hartree-Fock equations and find the coefficients that will take us to a molecular basis, we will see that the orbitals that are occupied in the atomic orbital basis, will have a contribution to occupied orbitals in the molecular orbital basis, and are therefore still important.

Contracted case

Lets now take a look at how the STO-$2$G basis set will be used in the case of contracted basis functions. Now our overlap matrix will take the following form:

$$ S= 
\begin{bmatrix}
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_1^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_2^\mathrm{H_2}\right.\right&gt; \\
\left&lt;\phi_2^\mathrm{H_2}\left|\phi_1^\mathrm{H_1}\right.\right&gt;  &amp; 
\left&lt;\phi_2^\mathrm{H_2}\left|\phi_2^\mathrm{H_2}\right.\right&gt; 
\end{bmatrix}
$$

At first glance the overlap matrix look just like the one in the uncontracted case, but smaller. What has changed is how the elements are defined:

$$\left&lt;\phi_i^\mathrm{H_k}\left|\phi_j^\mathrm{H_l}\right.\right&gt; = 
\sum_I\sum_J c_Ic_JN_IN_J\int_{-\infty}^\infty P_I\exp\left(-\alpha_Ir_I^2\right)P_J\exp\left(-\alpha_Jr_J^2\right) d\tau$$

We can see that we now sum over all the primitives included in the contraction. To make it more clear, let's write out how the second element in the overlap matrix looks like. Now remember that each of the basis functions consists of two primitives:

$$
\begin{eqnarray}
\left&lt;\phi_1^\mathrm{H_1}\left|\phi_2^\mathrm{H_2}\right.\right&gt;
&amp;=&amp; c_{11}c_{21}N_{11}N_{21}\int_{-\infty}^\infty \exp\left(-\alpha_{11}r_{11}^2\right)\exp\left(-\alpha_{21}r_{21}^2\right) d\tau \\
&amp;+&amp;c_{11}c_{22}N_{11}N_{22}\int_{-\infty}^\infty \exp\left(-\alpha_{11}r_{11}^2\right)\exp\left(-\alpha_{22}r_{22}^2\right) d\tau \\
&amp;+&amp;c_{12}c_{21}N_{12}N_{21}\int_{-\infty}^\infty \exp\left(-\alpha_{12}r_{12}^2\right)\exp\left(-\alpha_{21}r_{21}^2\right) d\tau \\
&amp;+&amp;c_{12}c_{22}N_{12}N_{22}\int_{-\infty}^\infty\exp\left(-\alpha_{12}r_{12}^2\right)\exp\left(-\alpha_{22}r_{22}^2\right) d\tau
\end{eqnarray}
$$

In this notation the first number refers to the basis function and the second refers to the primitive. As we can see each element of the overlap matrix in the contracted basis is made from four integrals. Even though our overlap matrix is four times as small compared to that of the uncontracted basis, we still have to evaluate the same number of primitive integrals. In this case we have investigated here, where we used STO-$2$G, we have no unoccupied orbitals. In general, this will not be the case, most basis sets (also when used in the contracted basis) will have unoccupied orbitals.

Other comments

In your question I can see that you mention molecular orbitals. Just to be clear. When we do the integrals we are in the atomic orbital basis. The molecular orbital basis is the one we find when we solve the Hartree-Fock equations. I.e. the molecular orbitals have nothing to do with uncontracted vs contrancted basis functions, though we will get more molecular orbitals if we use uncontracted basis functions, because we will get the same number of molecular orbitals as we had atomic orbitals.

Regarding how to implement this in Mathematica I have no idea. Though as my understanding of how Mathematica is ought to be used, I would recommend that you would implement your integral code in a language like Python. For inspiration you can take a look at this guide, which implements all the integrals with contracted basis functions.