A brief description of the techniques that we are using and general idea about my research can be found in this short review 

Andreev bound states in films (A. Vorontsov 2016)

The quasiclassical formulation of the Greens functions technique for many-body systems is one of the most powerful tools to investigate non-uniform states of superconducting condensates. The following review gives a good description of the technique with applications, is a great reference, but not so good as a first read. 

J. Serene and D. Rainer Physics Reports 101, 221-311 The quasiclassical approach to superfluid 3He

The best way to learn the technique is to start using it, and the simplest way to start is to self-consistently calculate the superconducting order parameter as a function of temperature. Just flipping through pages of the above papers will be enough to do that. Do the following steps for uniform superconductor:

  • Find the self-consistency equation on gap Delta (use Matsubara approach)
  • Find the quasiclassical propagator that enters the self-consistency equation
  • Eliminate cut-off and unknown interaction parameter from self-consistency
  • Solve the transcendental equation for Delta(T) numerically

For non-uniform superconducting states the simplest first project would be calculation of the gap profile near a specularly-reflecting boundary. Again, a simple search on the web will provide lots of references to compare your results. Take d-wave 2D material first. 

The self-consistency equation remains the same as before, except now we have to solve it on a spatial grid Delta(T,x)

The quasiclassical propagator has to be found by integrating Eilenberger Equation along straight trajectories with specular boundary conditions. It is better to not to attempt it directly, but rather to use parametrization of quasiclassical Green's function by coherence amplitudes. A detailed description of the `Riccati' parametrization can be found in 

Matthias Eschrig Phys. Rev. B 61, 9061 Distribution functions in nonequilibrium theory of superconductivity and Andreev spectroscopy in unconventional superconductors

This way one gets very stable Riccati-type differential equations for coherence amplitudes that are much nicer to integrate. The particular integration technique we are using, employs properties of Riccati equations, and is described here. One can probably also use some standard integration packages, but they might be slower. After the Riccati/coherence amplitudes are found for each trajectory, we use them to find the quasiclassical propagator for that trajectory. 

Insert that propagator into the self-consistency equation and again solve for Delta(T,x). The profile of the gap at the interface is given through x-dependence. Finding multi-component solution to self-consistency equation might be difficult, but one can use different tricks to accelerate this process. The one we are using is described here which is the Anderson Mixing method, see also

V. Eyert Journal of Computational Physics 124, (1996) 271 A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences