I am trying to solve an assignment on solving the Bogoliubov de Gennes equations self-consistently in Matlab. BdG equations in 1-Dimension are as follows:-
$$\left(\begin{array}{cc} -\frac{\hbar^{2}}{2m}\frac{\delta^{2}}{\delta z^{2}}-\mu+V\left(z\right) & \triangle(z)\\ \triangle(z) & \frac{\hbar^{2}}{2m}\frac{\delta^{2}}{\delta z^{2}}+\mu-V(z) \end{array}\right)\left(\begin{array}{c} u_{n}(z)\\ v_{n}(z) \end{array}\right)= \epsilon_{n}\left(\begin{array}{c} u_{n}(z)\\ v_{n}(z) \end{array}\right)$$ along with the equations for gap function $\triangle(z)$ and number density $n(z)$. $$\triangle(z)=U\sum_{n}\left(1-2f_{n,}\right)u_{n}(z)v_{n}^{\star}(z)$$ and $$n(z)=2\sum_{n}|{u_{n}(z)}|^{2}f_{n}+|{v_{n}(z)}|^{2}\left(1-f_{n}\right).$$
For the case of solving the BdG equations in Fourier space in Matlab for the case of a periodic potential and periodic gap function (assumed), we can take $$u_{n}(z)=\sum_{k}\exp\left[ikz\right]U_{n,k}, $$ $$\triangle(z)=\sum_{K}\exp(iKz)T_{K},$$ and $$ V(z)=\sum_{K}\exp(iKz)P_{K} $$ where the sum is over the reciprocal lattice vectors $K$ leading to the number equation $$N=2\sum_{n,k}\left[f_{n}|{U_{n,k}}|^{2}+\left(1-f_{n}\right)|{V_{n,k}}|^{2}\right]$$ with $ f_{n}$ as the Fermi distribution function. Solving the set of equations self-consistently for a fixed $N$, I am trying to get a value of chemical potential from the number equation each time after solving the eigenvector components $U_{n,k}$ and $V_{n,k}$, but due to the form of the exponentials in the number equation and sum over large number of them, I am unable to get a correct value of chemical potential out of them using Matlab routines as the root of the equation to put it back into the equations for eigenvector components.
In most cases, I get random values of chemical potential since the equation is more or less insoluble. How can I avoid this error ? Is there a better way to numerically solve the BdG equations self-consistently ? I also want to do this assignment in real space avoiding finite size effects but started with the Fourier space case to avoid errors associated with discretizing the differential. Please guide and ask for any details you might need.
Following is my MATLAB code to solve the equations in real space but the code does not work as fsolve does not find the mu value.
http://postimg.org/image/v7amx5vd9/full/
http://postimg.org/image/do9pyyk7z/full/
If you have solved BdG equations numerically, please tell about the method and steps you used such that the above problems are eliminated.
