With Mathematica it easy to do.

In Maple it's not so easy:

1 page.

2.page.

3.page

@ligonberry 

Maple have No  built-in collocation methods, least squares method, Galerkin method to solve ode's.

dsolve it have method `series`,but in your case dosen't work.

@radaar 

I made a mistake,Yes you are right for Z=4 sum is converge.

I check for i=100 and q=150.

PROBLEM1.mw

This really isn't a question in the absence of a concrete cut-and-pastable example,

and yours Maple 18  is NOT up to date?

@digerdiga 

From Maple's help:

_d01ajc is for finite interval of integration; allows for badly behaved integrands;

uses adaptive Gauss 10-point and Kronrod 21-point rules.

AND:

"maxintervals=300000". Probably is the maximum number of subdivisions for a given interval.BTW for [0,1].

That's all I know how much is written in the Maple help.

@kuwait1 

On my Maple 2018.1 works fine.

On MMA 11.3:

@kuwait1 

It's normal behavior if you want a converged sum you must  increase again q.

I do not understand , why you tried approximate the sum?, after all, you already know the answer: 0.1404283142

No_problem.mw

Increase Digits to 20 and q to 120.

Probably it's a bug.I'm submit possible bug to SCR.

@spalinowy 

Put your question in https://math.stackexchange.com and there they will tell you what 
I am telling you that it is impossible(No practical use)

Formula 15 will not help you because it is Laplace's Transform of the integral, 
and you have something else in the question.

@vv 

I copy code,paste,execute,dosen't work?

If I convert to 1D  then works. 

I'm confused.

@shimaa sadk 

On Maple 2018.1 yours attached file works fine.Maple calculate  all worksheet about 30 seconds.

See result:

Download aquestion_(2).mw

If you have problem then  try this worksheet below:

aquestion_(1)_ver2.mw

@shimaa sadk 

For general integral with general constans :

int(lambda[1]*alpha^2*exp(lambda[1]*Z)/((exp(lambda[1]*Z)-1+alpha)^2*(exp(lambda[2]*Z)-1+alpha)), Z = 0 .. infinity)

Maple 2018.1 can't find symbolic solution.

r[1] := proc (lambda1, lambda2, alpha) options operator, arrow; lambda1*alpha^2*exp(lambda1*Z)/((exp(lambda1*Z)-1+alpha)^2*(exp(lambda2*Z)-1+alpha)) end proc;
R[1] := int(eval(diff(r[1](lambda1, lambda2, alpha), lambda1), [lambda1 = 3/10, lambda2 = 2*(1/10), alpha = 4]), Z = 0 .. infinity, numeric);


#R[1] := 1.077997357

See attached file for more information:

aquestion_(1).mw

@spalinowy 

I gave you only a comment,my point of view. Wait for the answer.

And so at the end you'll end up with numerical calculations.

You have a Transcendental equations that why Maple can't solve(Any CAS can't solve).Transcendental equations have no analytical solution.

Try numeric solver like fsolve, DirectSearch package,or something else.

See:https://en.wikipedia.org/wiki/Transcendental_equation

        https://www.mapleprimes.com/questions/88377-Direct-Search-Optimization-Package