I often prefer to switch to fixed integration bounds. Then I get it to work.

Download MP-223127.mws

I am not very familiar with the G function, so the following is a bit formal.

Maple's notation with 'empty parameters' occures if some cancel out. Conversely one
can recover by feeding them (generically and in valid ranges).

The task means MeijerG([[a], [b]],[[c], [d]],z) with b=c.

Now take the series presentation in z=0 from the FunctionAdvisor, [b=c, d=1, a=1/2],
and ignoring possible issues (radius of convergence?)

The summands are -GAMMA(1/2+c+k)/Pi* sin(Pi*k) /GAMMA(c+k) * z^(c+k)*(-1)^k.

This is zero (if c is not a negative half-integer).

> [seq(2*n+1, n= -4 .. 4)];

                   [-7, -5, -3, -1, 1, 3, 5, 7, 9]

> map('q -> q mod 3', %);

                     [2, 1, 0, 2, 1, 0, 2, 1, 0]

>

Passing to numerators it becomes a polynomial system (check solutions against denominators later), use exact coefficients.

One can eliminate beta.

Then one has to solve a polynomial in delta, degree = 85 and coefficients of magnitude ~ 10^150.

From that one gets the beta values as a rational function of delta (coefficients also of high magnitude).

I have not carried it out in detail.

You may use [Re(g2), Im(g2)]; plot(%);

Your integrand has a singularity log(E)/sqrt(E) in 0, but you can start in E=1. Note that J:=Int(K, E=1 .. t) is / can be seen as a specific indefinite integral ("they all differ just by a constant"), which is 0 in t=0. So you can do the following:

# S = int(..., E) in the worksheet
K:=IntegrationTools:-GetIntegrand(S);
#series(K, E=0, 2) assuming 0 < E;

Int(K, E=1 .. t);
plot(%, t = 1 .. 10000);

J := -(exp(-a*s)-1)*kw/((exp(-a*s)+1)*s^2);
convert(J, trig): convert(%, tanh):
simplify(%): convert(%, tanh);

You may export as RTF (Word can read that) or as HTML (and pick the images from the according directory)

"If the file is being opened for writing and does not exist (and is not already open), it is created, if possible; otherwise, fopen generates an error.  If the file does exist (and is not already open), it is truncated. Writing starts at the beginning of the file if WRITE was specified, and at the end of the existing file if APPEND was specified."

Then you can use iostatus(), for example.

In my ancient times of coding at least a check for existence was done (and available handles), though I do not expect Maple to provide all the stuff for a coding laguage.

Download 222149.mw

It is a convention what 0^0 should mean and Maple treats it as follows:

`0^0`: %= parse(%);

                               0^0 = 1
`0^r`: %= parse(%);

                               0^r = 0
`0^1`: %= parse(%);

                               0^1 = 0
`0^(-1)`: %= parse(%);

   Error, numeric exception: division by zero

Actually you expect your result to be 11 * `0^0` = 11.

You can however use your own powering:

pow:= (t,n) -> piecewise( t=0 and n=0, 1, t^n);

m:=2; # or your m=10
x:=0;

Sum(f(x,j)*(m-j+1), j=0..m); # or your original double sum
value(%);
subs(f=pow, %);
eval(%);

                          2
                        -----
                         \
                          )   f(0, j) (3 - j)
                         /
                        -----
                        j = 0

                   3 f(0, 0) + 2 f(0, 1) + f(0, 2)

                3 pow(0, 0) + 2 pow(0, 1) + pow(0, 2)

                                  3

Int(sin(phi)*sin(-arctan((p*cos(phi)+Z)/sqrt((p*sin(phi))^2+r^2)))*hypergeom([5/2, 3/2], [5/2], sin(-arctan((p*cos(phi)+Z)/sqrt((p*sin(phi))^2+r^2)))^2)/((p*sin(phi))^2+r)^2, p = -1 .. 1);
tmp:=subs(phi=Pi/6, %); # = your function, evaluated
plot(op(tmp)); # "visually" correct

simplify(tmp);
   -32*3^(1/2)*(Int(p*(p^2+25)/((p^2+100)^(3/2)*(p^2+20)^2), p = -1 .. 1))

Now just argue by symmetry in zero.

Use fsolve( ..., complex)

As far as I can see you have a set of 36 quadratic equations. It is very likely that there is more than 1 solution (and it is not clear at all why those [complex ?] solutions are discrete). It is also possible that there is NO solution (avoiding "evalf" seems to indicate that).

PS: you can avoid the lengthy display by just uploading you sheet or a pdf of it.

rsolve({a(n) = 2*a(n-1) + 4*a(n-2), a(0)=1, a(1)=2}, a(k));

Your integrand has a singularity in 0.161144379880543e-3, so you can not use numerical integration directly. However the integral exists as Cauchy Principal Value (if it is what you want):

(E*sigma*(Int((-xi^2+1)*(-KImax^2+Kfc^2)/((KImax^2-Ksc^2)*((-R^4+1)*(KImax^2+Ksc^2)+eta*E*sigma)),
  l = l0 .. b, CPV))/alpha0);

value(%);

    8491.79571174594

where b is your upper bound 12.1*10^(-3)