This is the Fermi function evaluated at Tc

fermi[x_]  :=  1/(Exp[x]  +  1)

This is the BCS occupation probability vk^2, versus reduced energy x = ζ/(kBTc). ζ=ϵ-μ.  It assumes N(0)V = 0.42 (this gives a value for the cutoff energy hbaromegac of 9.00 kBTc).  We also take Δ(0) to be 1.76 kBTc (redgap = 1.76).  The plot goes from minus hbaromegac to plus hbaromegac.

vksq[x_,  
        redgap_]  :=  (1/2)*(1 - x/Sqrt[x^2 + redgap^2])

The Fermi function at Tc is virtually identical to the BCS occupation distribution at T=0 !!

Plot[ {fermi[x],  vksq[x, 1.76],  1 - vksq[x, 1.76]},  {x, -9, 
        9},  PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1],  
          RGBColor[0, 1, 0]}]

This is the difference between the Fermi function at Tc and the BCS occupation distribution.  They are the same to within about 2%

Plot[ fermi[x]  -  vksq[x, 1.76],  {x, -9, 9}]

Now look at the product uk * vk that appears in the definition of the gap function.

Plot[ 
      Sqrt[vksq[x, 1.76]]*Sqrt[1 - vksq[x, 1.76]],  {x, -9, 
        9}]

Plot[ 
      vksq[x, 1.76]*(1 - vksq[x, 1.76]),  {x, -9, 9}]