I have simplified my problem here, in my actual problem the matrix is much bigger which makes it impossible to find eigenvalues analytically. So, I chose standard BCS problem (2x2 matrix) to demonstrate my problem with Mathematica. But, with even this simple case, it takes forever to integrate. My code is the following:
hamiltonian[k_, mu_, delta_] := {{k^2 - mu, delta}, {delta, -k^2 + mu}}
eigens[k_, mu_, delta_] := Eigenvalues[hamiltonian[k, mu, delta]]
fermitotal[beta_, k_, mu_, delta_] := Block[{ee = eigens[k, mu, delta]}, 1/(1 + Exp[beta ee[[1]]]) + 1/(1 + Exp[beta ee[[2]]])]
nTotal[beta_, mu_, delta_] := NIntegrate[fermitotal[beta, k, mu, delta] k^2, {k, 0, 20}, Method -> {"LocalAdaptive", "SymbolicProcessing" -> False}, AccuracyGoal -> 4, PrecisionGoal -> 4, MinRecursion -> 10, MaxRecursion -> 300, WorkingPrecision -> 13]
nTotal[50,1,1/10]//AbsoluteTiming
Any idea how to speed up the calculation?
I have tried, using Compile for the hamiltonian matrix. Specifying the precision for the numbers,
nTotal[N[50,30],N[1,30],N[1/10,30]]//AbsoluteTiming
Did not help. I am completely stuck here and I don't know what causes this slowness. Any help would be very much appreciated.