@Preben Alsholm 

Thank you very much. It is the first time when I see the utilizing of the functiuon dsolve for the computation of the integral. Nice approach!

@Mariusz Iwaniuk 

Thank you for your help. I think the main problem is that Maple does not know the L'Hôpital's rule and how to solve the limit for the uncertainty 0*infinity :)

@acer 

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. 

restart;

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?

@vv 

Thanks for the help. I try to change infinities to the range n=-1000..1000 and it works. It gives me good approximation.

The problem is solved.

@tomleslie 

Sorry! I missed some part of the code. Here is the correct one

restart;
R0 := exp((2*Pi*I)*n^2*z);

z:=exp(I*Pi/k);

R1 := sum(R0, n = -infinity .. infinity);

R1 := abs(R1)^2;
R2 := exp((2*Pi*I)*(n+1/2)^2*z);
R2 := evalf(sum(R2, n = -infinity .. infinity));
R2 := abs(R2)^2;
R := evalf(sqrt(Im(z))*(R1+R2));
plot(R, k = 1 .. 10);

@_Maxim_ 

A lot of thanks. I tested the code. It works great!.

I apologize but may I ask you again. I try to extend your code for three variables and plot the solutions. But at the end I get error

@Carl Love 

Thanks for your response. Nevertheless Maple can solve this system of Eqs. for given value of t. See below

@Carl Love 

restart;

W1:=0.08687741457; W2:=2.584713564 ;W12:=0.5550746999; W21:=0.4045459677; 

R0 := x*W1-W12*y+ln(t)*x-evalf(Sum(2*Pi*t*(x/sqrt(((2*n+1)*Pi*t)^2+x^2)-x/((2*n+1)*Pi*t)), n = 0 .. 5000)) = 0:

R1 := -x*W21+W2*y+ln(t)*y-evalf(Sum(2*Pi*t*(y/sqrt(((2*n+1)*Pi*t)^2+y^2)-y/((2*n+1)*Pi*t)), n = 0 .. 5000)) = 0:

@ecterrab 

I forgot to say thank you. Sure, you are right. It was my stupid misunderstanding :)

@phil2  The first machine was the desktop with AMD FX 4100 quad-core processor (3.6 GHz). The second machine was based on AMD (as I mentioned before) but unfortunately I didn't fix the speed or the number of cores.

@phil2 Thanks a lot for your comment. Of course every time I used the 64-bit version of Maple. During last month I revealed very interesting empirical result related to this problem. I found out that there is no problem with numerical calculations if the machine is based on Intel proccessor instead of AMD one. I'm not sure that this is correct but I tested on 6 machines: 2 AMD-based and 4 Intel-based with 64-bit version of Maple.

@Markiyan Hirnyk 7178 Thanks a lot for the response, but I know that for some parameters this equation has no solution and anyway there was no problem to obtain the matrix on another desktop without lost of connection with kernel.

@Carl Love Thanks for the response and your answer. 

Let's take a look at the implicit function y=y(x) for different values of the parameter a=1..25

Tc_vs_thickness_beta=0_tau=9-975_25steps.pdf

We can see that there is an interval between approximately x=1 and 2 and x >3 (for some curves) where the function is not defined (i.e. piecewise defined). 

@Carl Love It's for my scientific calculations related to the physics of superconductivity. This code and transdental equation allows me to find the behaviour of the system, namely, critical temperature (denoted by y) in a dependence of the thickness x and some internal characteristics $\tau$ and a. On the next step I should understand via matrix, where this equation has no real roots (NaN) and then, of course, as you mentioned earlier I will plot the temperature distribution as a function of the thickness and the parameter a (contour plot). 

Btw since parameter a has a range from 1 to 50 this string in your code should look as

(i,j)-> R0f(i, 0 + (j-1)*(5-0)/(1000-1)),

instead of

(i,j)-> R0f(1 + (i-1)*(1-0)/(50-1), 0 + (j-1)*(5-0)/(1000-1)),

Is it correct?

@Carl Love It was exactly what I needed. Thank you very much!

@John May Also thank you for your attention to my question.