In the paper by Caroli, de-Gennes and Matricon on Bound fermion state in superconductor, they have reached a system of first order differential equation which reads roughly as: \begin{align} -i\hbar v_F\frac{df}{dr}+\Delta(r)g(r)&=(\epsilon+\frac{\mu\hbar^2}{2mr^2})f(r) \\ i\hbar v_F\frac{dg}{dr}+\Delta(r)f(r)&=(\epsilon+\frac{\mu\hbar^2}{2mr^2})g(r) \\[.5em] v_F&=\frac{\hbar k_F}{m} \end{align} $f(r),g(r), \Delta(r)$ are functions of $r$.
$\mu$ is an integer and $\epsilon$ is a small number and can be treated as a small perturbation.
What I fail to understand is how they are able to solve the above system by treating the right hand side as a small perturbation, to first order, and obtain the following results: \begin{align} f(r)&=e^{\frac{1}{2}i\psi(r)-K(r)}\\ g(r)&=-ie^{-\frac{1}{2}i\psi(r)-K(r)} \end{align} where $$K(r)=\frac{\int_{0}^{r}\Delta(r')dr'}{\hbar v_F}$$ and $$\psi(r)=-e^{2K(r)}\int_{r}^{\infty}e^{-K(r')}(\frac{2\epsilon}{\hbar v_F}+\frac{\mu}{k_Fr'^2})dr'$$ If someone can explain briefly the procedure to get to the final results, it would be very helpful.
Thanks in advance for anyone willing to answer.