seems to be too high ...

In your case however it is easy:

((x^7+1)^4/x^30); expand(%);

                     1     4     6     4      1
                    --- + --- + --- + ---- + ----
                     30    23    16     9      2
                    x     x     x      x      x

from which you can read it off or continue by

[op(%)];
map('t -> residue(t, x=0)', %);
convert(%, `+`);

                                  0

I do not know an actual answer to your question and I use M12.

a) Maple can be convinced to give a 'general' anwer

  Int(1/(x^4+3*x^2+1)^(p),x=0..infinity);
  Change(%, x=1/2*(-6+4*t)^(1/2),t); 
  # motivated by the symmetry of the *real* partial fraction decomposition
  # x^2+3/2=t, certainyl invertible on Positives
  value(%); combine(%);
  R:=%;

I omit to print the result, but Maple seems not to use reflection formula
to go further (I think it should be possible)

b) Now with Digits:=14 compare

  pTst:=1 + 1e-10;
  'eval(R, p=pTst)'; evalf(%);
  'Int(1/(x^4+3*x^2+1)^(p),x=0..infinity)': 'eval(%, p=pTst)'; evalf(%);
 
                           0.70258767772836

                           0.70248147304481

Is it really only correct for 3 leading decimals?

c) What happens here? Certainly a bug ...

  convert(R, StandardFunctions) assuming p::integer;

                                  0

I do not understand this today ...

It also may depend on your 'expression'. If it is not to complicated then it is better to post it as well.

You do not provide the situation, where the error occurs (the linked sheet works for me).

From the error message I would guess: you have negative values for x.

In your other post it seems that it can not be solved for beta ...

Let me write it in a more readable form, using x as the variable:

sqrt(Pi)*beta=-1/h*SteL/(exp(h^2*beta^2)*erf(h*beta))+
  R1*SteV/(exp((beta/R1)^2)*erfc(beta/R1));
subs(SteL=a,SteV=b, R1 = r, beta=x, %): 
B0:=combine(%, exp);
B0/sqrt(Pi)/x; 
B1:=collect(%, exp);

            1 = (some complicated stuff in exp and erf)

theAssumption:=a<0, b<0, h>0, r>0, x>0;

Now staring at it one can guess the pattern:

1 = '-b/2*D(erfc)(u)/erfc(u)/u - a/2*D(erf)(v)/erf(v)/v'; 
subs(v=h*x, u=x/r, %);
% = (B1);
is(%);
                                true

You can put the factor 1/2 into your a and b to have it nicer,
but that equation does not cry for an explicite solution in x.
I do not even see a series approximation.

So my feeling says: you need numerical solutions. 

Computing that will need some care: in x=0 you have a pole,

  B1; rhs(%);
  series(%, x=0, 5); because of D(erf)/erf

    (-1/2*a/h^2)*x^(-2)+b/Pi^(1/2)*r*x^(-1)+(1/3*a+2*b/Pi)+O(x)

and for x ---> oo that converges to b:

  B1; rhs(%); 
  eval(%, erf = erf(infinity)); # it gives 1
  asympt(%, x, 4) assuming 0 < r, 0 < h;
  lprint(%);

                     b+1/2*r^2*b/x^2+O(1/x^4)

but I have not looked inbetween (you should: if it is decreasing,
then there can not be a solution, if 1 < b=SteV ...).

But as alraedy said: I doubt there is a symbolic solution and it
is not even clear, whether there always exists a (real) one.

A:=Pi^(1/2)*beta = 1.878297101/exp(20.00000001*beta^2)/erf(4.472135956*beta)-
    2.100000000*R1/exp(beta^2/R1^2)/erfc(beta/R1);
identify(%): convert(%, rational):
A:=lhs(%) - rhs(%);

plots[implicitplot](A, R1=0..1, beta=0..1);


Now determine the 'extreme' points you are interested in (guessing from the graphic):

0=limit(A, R1=0, right) assuming 0 < beta;
expand(%);
beta1:=fsolve(%, beta=0.15);

                        beta1 := 0.1978948870

eval(A, R1=1);
beta0:=fsolve(%, beta=0.15);

                        beta0 := 0.1184172673


Then solve the equation 'formally' for R1 (not for beta, that gives troubles)
and consider that as a function of beta (should work due to the graphic):

RootOf(A, R1);
g:=unapply(%, beta);

Your graphic becomes much more smooth (it is the reflected one, correct that),
thus the original graphic would give quite error prone results:

plot(Re(evalf('g'(t))), t=beta0 .. beta1, numpoints=5);
plottools[reflect](%, [[0,0],[1,1]]);


If you want values then for example you can use

myStepSize:=0.01;
[seq( [t, g(t)], t=beta0 .. beta1, myStepSize)]:
evalf(%); 
                          myStepSize := 0.01


  [[0.1184172673, 1.000000000], [0.1284172673, 0.8441382672],

        [0.1384172673, 0.7054733505], [0.1484172673, 0.5804226649],

        [0.1584172673, 0.4661017791], [0.1684172673, 0.3598910320],

        [0.1784172673, 0.2587389755], [0.1884172673, 0.1567487861]]

You can find a recipe here: www.math.rwth-aachen.de/mapleAnswers/html/313.html

PS: it is worth to look at this pages beyond that question (even if it a bit old, some stuff may be outdated)

you can look at that in the online help (at the top of this page you find a search, activate 'help' and enter "DiscreteTransforms[FourierTransform]"

www.maplesoft.com/support/help/AddOns/view.aspx

The following might do, if you work with 'Standard' interface

plot3d(-x^3-y^2, x=-1..1, y=-1..1, axes=boxed, style=point);

If you work with the classical interface, then the graphics are less fine,
points can not be 'sized' (at least for me), but other symbol should do:

plot3d(-x^3-y^2, x=-1..1, y=-1..1, axes=boxed, style=point, 
  symbol=circle,  symbolsize=15, numpoints=20^2);

After simplifications this is a very special case, which Maple does not answer:

Sum((k+j)^j*r^j/j!,j = 0 .. infinity); eval(%,k=1);

                         infinity
                          -----          j  j
                           \      (j + 1)  r
                            )     -----------
                           /          j!
                          -----
                          j = 0

I googled for the keyword, seems somewhat special, both in combinatorics and other fields (epidemiology, queing ...)
and now think, that you have to invest some theory for that topic: search for "arch98v127.pdf" gives me the impression,
that your formula stems from a Lagrange inversion. And that cries to be treated before using it.

a) why?

b) what do all elements have in common and what can you operationally say for that?

There is no minimum. Even if I suppose it is homework:

The last terms are linear, where x,y,z live in a cube, thus this gives
something bounded (IIRC Analysis correctly).

So your 1st term matters. Transform to positive x,y,z:

  subs(z=-z,%);  subs(y=y+1, %);

    a*b^x*c^(y+1)*d^(-z)-x-2*y-2+3*z

now with 0 < z < 1, 0 < y < 2

This suggests you want b,c,d to be positive to avoid imaginary values
(may be you forgot to say that).

Now choose: subs(x=1/2, y=1, z=1/2, %);

This gives a*somePositiveConstant - anotherConstant

So for a ---> - infinity ...

Or similar ...

I do not have it in Maple, but I have snipped of the part up to the
'milestone 10', which gives the radix and U1, U2 for separating 1
(and to get the largest possible mantissa, that's what I wanted).

Not what you hoped, but probably small enough to get it in Maple without
too much pain (guessing that the Basic version would translate better,
but I do not have such). May be it helps a bit.

At least for evalhf I am sure, that it does not comply with IEEE 754.
As far as simulating compiled versions it would be nice to have such
a variant.

Download 102_separation.zip

a) You can use Wikipedia (or Google) to look up recursion or recursive

b) Hint: you want to look at something like this:

  The methods halfs an interval and depending on a decision you
  apply the method again for one of the halfs (until you stop).

gen:=proc(x0,a,b,M,n)
...
end proc; # this is missing in your code

BTW: if you use the Classical interface you may find to have less troubles starting with Maple,
for example it would have said "Warning, premature end of input", which tells you much more.