# Self-consistent numerical solution of a set of equations

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 $f_{n}=\frac{1}{\exp(\beta(\epsilon-\mu))+1}$. 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/

• It seems that this question is now gone from Physics as well. So you can leave it here. Do keep in mind though that cross-posting within the network is something that is frowned upon and usually not tolerated. Certainly not by posting just the link, as you did before. Just a friendly warning not to do that, because getting suspended because of it would be annoying. Just keep your questions limited to a single site, even if it fits on multiple ones. Good luck. – Bart Mar 9 '14 at 11:16
• Hi @user38579 and welcome to scicomp! Could you elaborate as to why you suspect "the form of the exponentials" is a problem? How many terms do you have to sum up, typically? How large are these values that you're adding up? – Paul Mar 10 '14 at 14:48
• @Paul : thanks. I have to typically sum over a 1000 distribution functions containing exponentials. Fourier Coefficients in the sum may have values as the fourier coefficients of $sin^2(z)$ with 2m/(h*a)^2 of the order of 10^30. Now that you asked this question and thinking about how froot works, I am thinking that problem is not with the exponentials but rather might be with the size of coefficients- but I thought Matlab could handle 10^30.Instead of multiplying 2m/((ha)^2) as I did to work with large numbers >1 , should I multiply by 2m/((ha)^2)*(10^15) ? I hope you get what I am saying. – user38579 Mar 11 '14 at 15:02
• Is that $\Delta(z)$ in the lower left supposed to be transposed? – Nick Alger Mar 17 '14 at 14:09
• @NickAlger : it is a diagonal block matrix when written in matrix form since it has no derivative term. In general, it is conjugate of the upper-right but we can assume $\triangle(z)$ to be real and thus lower-left and upper-right are same. – user38579 Mar 17 '14 at 20:27