Differentiating a function that has NIntegrate within it

Question

I am something of a novice with Mathematica. I have written the following code;

Clear[band, en, w, fermi, k, T, S, a]
k = 86*10^-6;
T = 4000;
band[en_, w_] := 10000 Exp[-en^2/2 w^2];
fermi[en_, ef_] := 1/(Exp[(en - ef)/(k T)] + 1);
Plot[band[en, w] /. w -> 1, {en, -4, 4}]
Plot[fermi[en, ef] /. ef -> 0, {en, -4, 4}]
Plot[{band[en, w] /. w -> 1, (band[en, w]*fermi[en, ef]) /. {ef -> 0, w -> 1}}, {en,-4,4}]
S[ef_] := -k NIntegrate[(band[en, w]*fermi[en, ef]*Log[fermi[en, ef]]) /. w -> 1, {en, -Infinity, Infinity}]
Plot[S[ef], {ef, -4, 4}, PlotRange -> All]
a[ef_] := D[S[ef], ef]
Plot[a[ef], {ef, -4, 4}]

It works very nicely until I try to plot a, which is simply the derivative of S wrt ef - even though I can plot S vs ef a few lines above. Plotting a gives error messages along the lines of;

General::ivar: "-3.83657 is not a valid variable."

This is clearly something to do with NIntegrate evaluating too early, or not being a function of ef in the way I would expect it to be, or similar. But I have no idea how to fix it.

I would be very grateful for any suggestions anyone could give. I think this issue is showing up my fundamental lack of understanding how Mathematica works, and so replies along the lines of "You should read and understand such-and-such article..." would also be very welcome!

Many thanks,

Nick

Answer 1

Use a numerical derivative:

Clear[band, en, w, fermi, k, T, S, a]
<< NumericalCalculus`
k = 86*10^-6;
T = 4000;
band[en_, w_] := 10000 Exp[-en^2/2 w^2];
fermi[en_, ef_] := 1/(Exp[(en - ef)/(k T)] + 1);
S[ef_?NumericQ] := -k NIntegrate[(band[en, w] fermi[en, ef] Log[fermi[en, ef]]) /. w -> 1, 
                   {en, -Infinity, Infinity}]

Plot[{S[x], ND[S[ef], ef, x]}, {x, -4, 4}]