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!

Kevin