A  and  B  are 2 curves on the surface  C :

restart; 
f:=(x,y)->2*x^2+3^2-x*y-4:
A:=plots:-spacecurve([1,y,f(1,y)], y=-2..2, color=red, thickness=3):
B:=plots:-spacecurve([x,1,f(x,1)], x=-2..2, color=blue, thickness=3):
C:=plot3d(f, -2..2, -2..2, style=surface):
plots:-display(A,B,C, orientation=[115,75], labels=[x,y,z]);

For this particular plotting, 20 digits is enough. It is better not to specify the global Digits, but to specify the desired accuracy when specifying the function itself. In this case, you need a procedural form for plotting:

restart;
f:= x->evalf[20](exp(x)*HeunC(2*x, 1/2, -1/2, -x^2, 1/8, 99/100)):
plot(f, 0..40);
                          

restart; 
a := plot(1/x, x=0..1);
op([3,2], a);   # This is the desired result
 

Should be  od  (instead of  end )  in the end of your code.

You wrote "does computing H(n) for n greater than 10^16 give  any information about the series" . The convergence or divergence of the series is related to the behavior of this series at n->infinity. Therefore, without this study, computations with any accuracy of a finite number of terms of the series do not give an exact answer to the question of its convergence.

See help on  size  option for  plot  command.

Example:
plot(sin(x), x=-Pi..Pi, size=[800,300]);
plot(cos(x), x=-Pi..Pi, size=[800,300]);

If you need to enumerate in this long list in a for-loop, you do not have to first generate the list itself. You can do this dynamically using the following commands: combinat[nextperm], combinat[prevperm], combinat[firstperm] and combinat[lastperm]
 

Here is a simple way (it works in Maple 2018.1):

restart;
with(Grading):
Quiz("What is the sum", proc(Resp, Ans) evalb(value(Ans) = Resp) end proc, proc() local r1,r2; r1 := rand(1 .. 20)(); r2 := rand(1 .. 20)(); r1%+r2 end proc);
 

restart;
Eq1:=x+2*y+3*z=9:
Eq2:=4*x+y+5*z=29:
Eq3:=2*x+2*y+5*z=21:
solve([Eq1,Eq2,Eq3], [x,y,z]);
eval(2*x+3*y+z, %[ ]);

Below is another way:

restart;
Eq1:=x+2*y+3*z=9:
Eq2:=4*x+y+5*z=29:
Eq3:=2*x+2*y+5*z=21:
solve([Eq1,Eq2,Eq3, A=2*x+3*y+z], [A,x,y,z]);
 

Another method that is commonly used for manual calculation is to multiply and divide your fraction by the sum of the roots in the numerator:

Download rationalize_new.mw

Here is another solution:

restart;
N:=100:
L:=NumberTheory:-Divisors(N);
add(a(k), k=L);
                                 L := {1, 2, 4, 5, 10, 20, 25, 50, 100}
         a(1) + a(2) + a(4) + a(5) + a(10) + a(20) + a(25) + a(50) + a(100)


I think for large  N  this will be more effective than the method above.

1. If we expand  cosh(x)  in the power series, then we see that all the coefficients of the series are positive numbers and  x  has even exponents:

series(cosh(x), x=0, 10);
     1+(1/2)*x^2+(1/24)*x^4+(1/720)*x^6+(1/40320)*x^8+O(x^10)

Therefore, it is easy to obtain an entire chain of inequalities  (in real domain):
cosh(x) >= 1+(1/2)*x^2 ,   cosh(x) >= 1+(1/2)*x^2+(1/24)*x^4   and so on.

2. Since for any real  x  this is obvious that  abs(cosh(x))=cosh(x) , all the above inequalities are also hold for  abs(cosh(x)) , for example  abs(cosh(x)) <= exp(x^2/2)  and so on.

I guess that you are looking for a solution in a real domain. Since your equation has the parameter  k , then to get the full solution, use the options  parametric  and  real . The desired result (the positive root) can be obtained as follows:

solve(7/(x^2+1)=k, x, parametric, real);
A:=op(4, %);
B:=eval(x, select(t->sign(op([1,2], t))=1, A)[]);
sqrt(simplify(B^2));
              

restart: 
int( exp((a+b*I)*exp(I*x)) , x=-Pi..Pi)  assuming real;  # a  and  b  are assumed to be real
simplify(%)  assuming real; 
                                      2*Pi*signum(a)^2
                                                 2*Pi 

restart;
A:=2*a+2*b;
content(A)*``(primpart(A));
expand(%);
                           

See help on these commands for details.