you can continue by adding the following command

  subsindets(%, 'integer', evalf);

and if you dislike floats as exponents then moreover

  identify(%);

Note however, that this still stands as a compact form for
a polynomial to be solved (and since the numbers are
huge it may be a good idea to increase Digit).

Int = shows you the integral

int = evaluates it

I prefer to use

  Int(...)
  value(%)

Try it with sinus, for example.

The same with Diff and diff, Sum and sum, Product and product

Writing x=56/t and a simple re-scaling shows, that you are asking about solving
c = f(t) for the function f:= t -> arccosh(t)/t.

  Digits:=15;  
  f:= t -> arccosh(t)/t;
  plot(f(t), t= 1 .. 20);
  
  'limit(f(t),t=1, right)': '%'=%;
  'limit(f(t),t=infinity)': '%'=%;

That shows: the (continous) function usually will have 2 solutions as it approaches 0 in the
end points, except the (unique) extremum.

  '0=D(f)(t)'; 
  t[max]:=fsolve(%);
  'f(t[max])': '%'=%;
 
                      t[max] := 1.81017058069898
                    f(t[max]) = 0.662743419349184

For your task you want 40/56 = f(t), which is beyond the maximum, so there is
no real solution (as it was already said).

There is no 'known' symbolic (piecewise) inverse, so one has to work numerically.

If your constant is not already the maximum, then you can proceed as follows:

  c:=0.2; # just as an example
  fsolve(f(t)=c, t=1 .. t[max]), fsolve(f(t)=c, t=t[max] .. infinity);

                  1.02091799400673, 17.8803266646570

This even works for quite small c (i.e. large t), even if the function becomes 
very flat towards infinity.

But for large t one can do better:

Using convert(f(t), ln) shows that f(t) = ln(t + sqrt(t^2 - 1))/t and since the
'-1' numerically can almost be ignored in that region one arrives at 

  const=ln(2*t)/t;
  RootOf(%, t); allvalues(%); 

    -1/const*LambertW(_Z,-1/2*const)

Here _Z (or similar) stand for "some integer". Some experimentation show that
the correct solution towards infinity (and const towards zero) should be

  t = -LambertW(-1, -const/2) / const;

You may wish to check 0 ~ f(-1/const*LambertW(-1, -const/2)) - const by either
evaluating in test values or by plotting (but use enough Digits). 


Edited to add: for the other 'end', t ~ 1, one can solve const = f(cosh(u)) for
positive u without problems after simplification (since division by t = cosh(u) 
is more or less equal one). That avoids the problem with D(f)(1+) = infinity,
just plot to see what is meant.

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.