@vv  Thank you for the clarification. So the correct approximation of EllipticF(z,k) for large values of z is -EllipticCK(k)*I, right?

@Preben Alsholm a lot of thanks for your comment and your efforts. Frankly speaking I do not believe that a solution exists but I tried to find the interval of the parameter q where I can find roots. I worried that I do something wrong. That's why I decided to ask the community. 

@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).