with(inttrans);
h:= 15/(s^3+6*s^2+15*s+15);

invlaplace(h,s,t);

           -----
            \                                2
             )      (10 + 7 _alpha + 2 _alpha ) exp(_alpha t)
            /
           -----
        _alpha = %1

                                 3       2
                  %1 := RootOf(_Z  + 6 _Z  + 15 _Z + 15)

allvalues(%): evala(%);

That wil give you a lengthy expression, which you also
can see numerical, use evalf[n] instead.

Maple decomposes h as a partial fraction and more or
less computes invlaplace of factor * 1/(a-s)

You can try fsolve(equation, a = 0 ... infinity) and setting Digits high enough before.

Usually that works.

You can try evalc(Re( expr ) ), same with imaginary part.

I do not have maple at hand now, but may be you want to re-write
the inert Diff as D(p)(T) + I * D(q)(T), if the suggestion does not work

Fa:=(Fa1*cos(wt))*cos(theta);
Fb:=(Fa1*cos(wt-120*(2*Pi)/360))*cos(theta-120*(2*Pi)/360);
Fc:=(Fa1*cos(wt-240*(2*Pi)/360))*cos(theta-240*(2*Pi)/360);

Ftotal:=Fa+Fb+Fc;

expand(%): combine(%);

                       3/2 Fa1 cos(wt - theta)

The coefficient at the 1st and the 2nd derivative is the same and after dividing
the equation has the following form (abreviating your Phi[4][1] as Phi)

  DE := diff(Phi(r),r$2) + diff(Phi(r),r)/r -a*Phi(r)/r^2 -b*Phi(r) -c = 0;

Note that a = 36, but that depends on the working precision.

Maple has an answer for that through dsolve

  BesselJ(a^(1/2),(-b)^(1/2)*r)*_C2 +
  BesselY(a^(1/2),(-b)^(1/2)*r)*_C1 -
  c*LommelS1(1,a^(1/2),(-b)^(1/2)*r)/b

For a=36 the result would be

  BesselJ(6,(-b)^(1/2)*r)*_C2 + BesselY(6,(-b)^(1/2)*r)*_C1-
  c*(-23040+1152*b*r^2-36*b^2*r^4+r^6*b^3)/b^4/r^6

In former times people where able to ask questions in complete sentences.

And explained what they want (while everybody accepted that English is not
everybodys language).

Most of them also provided code for copy & paste & run.

It seems that times changed to fire-and-demand-reply.

How sad. And not the only case ...

I do not think, that it is in the set of characters for 'fprint'
(roughly: ASCII, like printing from a C program).

Without Export I get an error:

Error, (in LinearAlgebra:-LA_Main:-Eigenvectors) cannot determine if this expression is true or false:
abs(Re(h2))+abs(Im(h2)) < -.1800000000e-1+1/10*abs(Re(h2))+1/10*abs(Im(h2))

What is h2 ?


Besides that I would try first to look at results and use only some looping steps.

May be I even would write that to temporary space before exporting, i.e. before crossing system borders.

http://www.math.ufl.edu/~fgarvan/papers/qmaple.pdf

http://www.math.ufl.edu/~fgarvan/qmaple/1.0/maple11/win/wmprog

Have you already followed the download page and the instructions ?

http://www.math.ufl.edu/~fgarvan/qmaple/download_1_0.html

He also provides a tutorial covering that, in the link you gave

Digits:= 15: # to use fast numerical integration methods
g:=Delta -> exp( -2*
    Int(1/((x^2+25/4*Delta^2)^(1/2)*(exp((x^2+25/4*Delta^2)^(1/2))+1)),x = 0 .. infinity,
      method = _d01amc) # see the help about int/evalf
    ):

plot(g, 0 .. 1);        # unfair - but now very efficient ...

plot('g'(x), x=0 .. 1); # works as well

One can try to determine the 'peak' for the summands, but that - for me - does not
work in a satisfactory way (not even for using log).

But here is a different suggestion: asymptotically Zeta(q) = 1 + 1/2^q + 1/3^q +...

Depending on desired numerical precision one can find q, such that Zeta equals 1
(numerically), for double precision that would be given by DBL_EPSILON, q=52,
which is for p = q+2 = 54

Now approximate as

  'Sum(Zeta(p-2)*a^p/GAMMA(1/2*p+1),p = 4 .. infinity) =
 
   Sum(Zeta(p-2)*a^p/GAMMA(1/2*p+1),p = 4 .. 54) + 
    sum(   1    *a^p/GAMMA(1/2*p+1),p = 55 .. infinity)';

The last expression has a symbolic solution (in exp, the error function and some
polynomial in a), the first one is a finite sum (may be add should be used).

So one proceed by 

  eval(%, a=8): collect(%, exp);
  evalf[30](rhs(%)); # ~ 0.12 * 10^29 

             12470298161623233765889590075.3

If one wants to take more care one uses

  eval(%, a=8): collect(%, exp):
  op( rhs(%) );  

and after evalf recognizes a possible problem of numerical cancellation (terms
of similar magnitude to be added, but of different signs).

That could be 'healed' by increasing the working precision (just try to see it)

PS: guess the approach can be refined to smaller p, as the result is large.

Additional to Preben's anwer I miss a reply, what the parameters are: (positive?) integers?

And: what is the reason/application for the integral?

If it is pure exercise then this is not much motivation ...

Your function c(k) writes without almost any modification as 

  F := k -> 
    Int(Int(2*(-3*x+7*t)^k/GAMMA(k+1)*exp(x-7*t),x = 0 .. t),t = 0 .. infinity); 

                 infinity    t
                /           /                k
               |           |   2 (-3 x + 7 t)  exp(x - 7 t)
    F := k ->  |           |   ---------------------------- dx dt
               |           |           GAMMA(k + 1)
              /           /
                0           0

This is done by switching to uppercase Int (= 'inert') and combining exponentials

  subs(int = Int, c(k)); combine(%, exp); F:=unapply(%, k);

Instead of using the somewhat more simple version H(k) already given (missing any
reply, test results or confirmation - so I guess statements are enough for you)
finally I arrived at a complete solution. At least I think so. And just check on
values of those k, for which you are able to compute, i.e. c(k) = F(k) = f(k) = H(k).

  f:= k -> 1/7*(1 - (2/3)^(k+1));

But then summing that over k can not converge.

Edited: the sheet attached MP_2dim_integral_erl.mws

I do not have Maple at hand right now, so can not check on input errors
(which may occur using the standrad interface)

But in the displayed code you set x:=0, while it also is a variable for integration.

That would explain the shown error:

you want x = 0 ... t but the system gets 0 = 0 ... t