Finally got things working pretty well.

Plot 16 was a nice result, generated by 
>> Master(10,9,-5,1,20,1,200); 
Between iterations 130-180 it shows a halo of the excess particles around a BCS condensate. For the other iterations, n_up==n_down is stuck. 
The getting stuck thing seems to happen a lot with strong interactions (large U). Unsticking them requires mu_down to get very small, near 1, which is strange considering that it should be up around 8 or 9 in the non-interacting case.

As we increase the temperature we find that the condensate starts to dissolve. The critical temperature Tc is about the same as U for small polarizations.

For polarizations like (10,9), U=-5, kT=3 condenses, kT=4 barely condenses (Delta increasing at iteration 200), and kT=5 does not condense (Delta -> 0).
For polarizations like (10,7), U=-5, kT=3 condenses, while kT=4 barely does not condense (Delta decreasing at iteration 200).  
For polarizations like (10,2), U=-5, kT=3 does not condense at all.

Higher polarisations thus make it harder for _any_ condensation to occur.  

Condensation either happens or doesn't happen - you can't get stable Delta=0.2.