This is actually exactly the same question I was asking myself a while ago, and it took me quite some time to figure it out. 
What I ended up doing was:

 - diagonalise the BdG Hamiltonian in Mathematica
 - solve the expressions for the eigenvalues for $\kappa$ 
 - neglect terms proportional to $\kappa^2$, $\kappa \Delta$ in the expressions for the eigenvectors, expand them to first order in $\kappa$, insert the expression for $\kappa$ you got from the eigenvalues
 - simplify everything (using the expression for $\epsilon$ and approximate $\hbar v k \approx E_F+U$) and arrive at beenakkers simple looking expressions (tedious)

I hope that helps some, in case you haven't figured it out yourself in the meantime. If you need more info/details don't hesitate to ask, I can elaborate some more if needed!