Near the core u looks like u0*(1+u0*x^2) with u0=1/R^2
and R=K(Tc-t)^(0.5+alpha) with alpha(t)=1/sqrt(-2*log(Tc-t)).
Here we plot alpha(t) and R(t) (for Tc=0.5):
Comparison between quad (green) and double (pink) precision runs.
ubar (approx inside the core) in red
psi (approx outside of the core) in blue
The largest curves (resp. double and quad) corresponds to highest
possible u0 obtained numerically.