@Rouben Rostamian  

The problem has a simple physical background. In the solid-state physics, it is possible to determine the real and imaginary parts of the dielectric permittivity \epsilon through the polarization function, depending on the frequency omega and a certain parameter q. In fact and roughly speaking, the integral is the same polarization operator, which is defined through the law of dispersion . 

Usually, if one knows the dispersion law (e.g., the most common example is quadratic one), one can predict the frequency behavior of the dielectric permittivity function. I want to solve the inverse problem.

Let us assume that we know the form of the dielectric permittivity. Let us imagine that it has the form of a slowly varying function (https://en.wikipedia.org/wiki/Slowly_varying_function). The task is to determine the corresponding dispersion law. 

I must confess that I chose the integration limits almost arbitrarily. In principle, there are no physical restrictions on other values of integration limits. I wanted to understand for a start how it is technically possible to numerically solve this equation, which can be called as "delay-integral equation". 

@Carl Love 

Sorry for this misunderstanding. I remember you from your helpful responses to my previous questions, so I just wanted to say thank you. I really didn't know that this type (DDE) of equation exists. 
And here I am just looking for how to numerically solve this equation, or at least get a couple of hints that will push me towards solving

@Carl Love 

As usual, thank you very much for your comment and clarification. I had assumed that this type of equation referred to integro-functional equations. 

Above, I provided some additional information about the values of the parameters and the unknown function.

@Rouben Rostamian  thanks for your response.

1. You are right absolutely.  At least, at the begininig I wanted to undertand qualitatve behaviour of this function .

2. I forgot to mention that definitely I was seeking the numerical solution. It would be great if the analytical solution exists even in the form of very cumbersome expression or even in the form of highly nonlinear differential equation. 

As for numerical values of parameters:

a=1, b=100

and  is positevely defined. 

Values of omega are also positive.

@Preben Alsholm 

Thank you fo the answer and for the suggestion. Here for this differential equation a parameter p is an ineteger number. At least for p=0 there is an analytical solution expressed in terms of periodic functions, namely sine and cosine fucntions. So obviously based on that periodic conditions I expect to obtain periodic solutuion also.

I have tried already to use initial function as sin(x)+cos(x) but now I experience another problem "Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging".

@nm Many thank for the reply and for the improvment! It's just my unwanted error with the incorrect call of the type for the solution.

I have chosen n(x) as a function for the convinience to export later this expression to latex. 

@dharr 

Thank you for the comment and for the help. I found the solution (see my comment above).

@Preben Alsholm 

It's ok. I have resolved that issue and I have extracted the data.

Thanks a lot for the reply.

I know that the execution of the code is rather long but my problem is that I can't extract the data file from the plot due to this warning.

@vv  I apologize for the question again where did you find the asymptotic expansion of incomplete elliptic integrals in the case complex argument and modulus? I mean in what textbook? 

@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?