Wow! A lot of thanks! It looks amazing.
May I ask you additional question. My goal is to plot more complicated function (see the code below), which is represented by twice integration and subsequent summation.
tt := -2; T_c := 0.169064e-1; mu := .869262; k := 2;
Omega := 2*Pi*N;
R0 := a*tanh((a^2-mu)/(2*T_c))*ln((2*a^2+2*a*q+q^2-2*mu-I*Omega)/(2*a^2-2*a*q+q^2-2*mu-I*Omega))/q-2;
R1 := int(R0, a = 0.1e-2 .. 100):
R2 := evalf(int(q*ln((-q^2-k^2+mu+I*(2*N*Pi-w)+k*q)/(-q^2-k^2+mu+I*(2*N*Pi-w)-k*q))/(k*(tt+evalf(R1))), q = 0.1e-2 .. 100));
R3 := evalf(Sum(R2, N = 0 .. 100)):
plot(Re(R3), w= 0.1e-2 .. 10);
As I understood my main problem is the integration procedure, namely the result of twice integration is not a number... Could you explain me, what was wrong?