It makes sense to use Maple only for numeric purposes here. 
For the two examples, a few numeric computations will suggest the general result, but the proof must be given by hand.

a) The limit  (M⊗N)(a,b) = sqrt(a*b). The proof is simple, as mmcdara showed.

M := (a,b) -> (a+b)/2:
N := (a,b) -> 2*a*b/(a+b):
MN:=proc(M::procedure, N::procedure, a, b, n::nonnegint, nmax::nonnegint:=n)
  local an:=a, bn:=b, ab:=Vector(1), k; 
  ab[1]:=[a,b]; 
  for k to nmax do  an, bn := evalf(M(an,bn)), evalf(N(an,bn));
    if k>=n then ab(k+1-n):=[an,bn] fi 
  od;
 seq(ab);
end:
MN(M,N,1,2, 2,5);             # the sequence [a1,b1],[a2,b2], ...

    [1.416666666, 1.411764706], [1.414215686, 1.414211438], [1.414213562, 1.414213562], [1.414213562, 1.414213562]

MN(M,N, 1,2, 10) = sqrt(2.);  # [a10,b10]

            [1.414213562, 1.414213562] = 1.414213562

b)  The limit M⊗N is the logarithmic mean; this can be proven similarly to the Gauss' proof for the AGM mean.

L := (a, b) -> (b-a)/ln(b/a):  # the logarithmic mean
a:=2.:  b:=8.:
MN( (a,b) -> (a+sqrt(a*b))/2, (a,b) -> (b+sqrt(a*b))/2,  a, b, 50 ) = L(a,b);

        [4.328085124, 4.328085124] = 4.328085123

The correct spherical parametrization is:

plots:-spacecurve([5^(1/2)*(t^2+1)*abs(t), arctan(4*t^2,2*t^3-2*t), arccos(1/5*signum(t)*5^(1/2))],
t = -1 .. 1, coords = spherical);

(Note the order [r, theta, phi] in plot3d)

It seems that the types `*`  and  `&*`  have problems for polynomial expressions.

restart;
type(x^2*y,  identical(x^2) &* name );
#                             false

type(x^2*y,  '`*`'({identical(x^2) , name}) );
#                             false

type(x^a*y,  identical(x^a) &* name );
#                              true

type(x^a*y,  '`*`'({identical(x^a) , name}) );
#                              true

1.0/3  is automatically simplified to 0.333...   with the current value of Digits (due to the so called Floating-Point Contagion).

Compare:

Digits:=3; evalf[8](1.0/3);   # Floating-Point Contagion
                          Digits := 3

                             0.333

Digits:=3; evalf[8](1/3);  # No contagion, 1/3 is exact
                          Digits := 3

                           0.33333333

The integral is computed by ScalarPotential  (if it exists, i.e. the form is exact).
Here is an example and a calling procedure (it works for more than 2 variables).

int1:=proc(F::list, X::list(name))
  uses VC=VectorCalculus;
  VC:-SetCoordinates('cartesian'[X[]]);
  VC:-ScalarPotential(VC:-VectorField(F))
end:

int1([2*x*y^7, 7*x^2*y^6], [x,y]);
      x^2*y^7

If you want inert integrals, it's easy to modify ScalarPotential.

It seems that the multiplication of float matrices with dimensions > 2000 crashes Maple 2021.

restart:
n:=2001:  # n:=2000 is ok
A:=Matrix(n, (i,j)->(i+2*j)/1000., datatype=float[8]):
A^2: #crash

Ortho.mw

A basis always exists, but the coefficients are probably huge and many primes are needed.

Symbolic computation with floats is inherently problematic.

f:=piecewise(x<0, exp(1), exp(1.)):
diff(f,x);

             piecewise(x = 0., Float(undefined), 0.)

restart:
a:=1:
before := {anames('user')};
#                         {a}
b:=2:
after := {anames('user')} minus eval(before,1) minus {'before'};
#                         {b}

You want

series(f/f1, a=0, 7);
expand(%);   # optional

but you have missing terms and other mistakes in your expected result. 

Your system is polynomial. I have substituted g = G*K.
The system reduces to find the roots of a huge polynomial in I1 of degree 9, with 8 parameters.
There will be 9 solutions. A few hundred pages will be enough for the result (if you arrange them nicely).

restart;
F:=[b*R + b*S - d*S - bcd*I1*S - bcd*I2*S - (e*g*S)/K - (g*I1*S)/K - ( g*I2*S)/K - (g*R*S)/K - (g*S^2)/K,
-d*e - (e^2*g)/K - (e*g*I1)/K - (e*g*I2)/K - (e*g*R)/K + bcd*I1*S +  bcd*I2*S - (e*g*S)/K - e*scd1,
-d*I1 - (e*g*I1)/K - (g*I1^2)/K - (g*I1*I2)/K - (g*I1*R)/K - ( g*I1*S)/K + e*scd1 - e*f*scd1 - I1*scd2,
-acd*I2 - d*I2 - (e*g*I2)/K - (g*I1*I2)/K - (g*I2^2)/K - ( g*I2*R)/K - (g*I2*S)/K + e*f*scd1,
-d*R - (e*g*R)/K - (g*I1*R)/K - (g*I2*R)/K - (g*R^2)/K - (g*R*S)/K +  I1*scd2]:
P:=eval(F, g = G*K):
indets(P);
#      {G, I1, I2, R, S, acd, b, bcd, d, e, f, scd1, scd2}
V:=[S,e,I1,I2,R];
#                     V := [S, e, I1, I2, R]
eliminate(P,{S,e,R}):
sys:=[%[2][]]:
map(degree, sys, [I1,I2]);
#                             [5, 5]
eliminate(sys,{I2}):
polyI1:=%[2][]:
degree(polyI1,I1);
#                               9
nops(polyI1);
#                               10
length(polyI1);
#                            507655
convert(polyI1,string): length(%);
#                             432373
indets(polyI1);
#             {G, I1, acd, b, bcd, d, f, scd1, scd2}

Maple is not able to recognize a particular 0 expression. This is not uncommon in any CAS. In our case:

ex := -sqrt((-8*x - 16)*exp(x/2) + x^2 + 4*x + 16*exp(x) + 4) + 4*exp(x/2) - x - 2:

is(simplify(ex)=0) assuming x::real;
        FAIL
simplify(eval(ex, x=2*t)) assuming t::real;  #workaround
        0

I don't think this is a true bug. The designer simply did not anticipate that a user would choose numerical labels. Probably such numerical labels should not be accepted.

You cannot prescribe all three points of tangency, because this implies 6 conditions and the ellipse has only 5 parameters.
But it is possible (in general) to fix two tangency points and the third will be computed.