I couldn't find a stock method for this, so I wrote a procedure. The procedure will work for any 2D PLOT structure containing two or more CURVES substructures.

ShadeBetween:= proc(
    P::specfunc(PLOT), 
    {pairs::list(patlist(posint, posint, {name, name= anything})):= [[1,2]]}
)
uses PT= plottools, LT= ListTools;
local 
    C:= subsindets(
        op~(1, select(type, [op](P), specfunc(CURVES))), rtable, convert, 'lislist'
    ),
    p
;
    plots:-display(
        P,
        seq(PT:-polygon([C[p[1]][], LT:-Reverse(C[p[2]])[]], p[3..][]), p= pairs),
        _rest
    )
end proc
:
#Example usage:
A := <1, 2, 3, 4, 2, 3, 1>:
B := <1, 5, 1, 0, 2, 1, 1>:
C := <1, 2, 1, 3, 1, 4, 2>:
ShadeBetween(
    Statistics:-LineChart([A, B, C], format= stacked), 
    pairs= [[1, 2, color= pink], [2, 3, color= orange, transparency= 0.3]]
);

You can't just substitute vectors for scalars in arithmetic expressions; however, Maple's elementwise operators (indicated by the ~ character) make it almost that simple. You can handle this very similarly to your previous Question:

restart:
N := 2.7:
Hb := <0.076, 0.083>: k := <0.00003400566801, 0.00003424620533>:
P50a := <20.78938475, 21.39546041>:
P50v := <21.20711722, 22.06793197>:
nu := <0.02042461957, 0.02120393111>:
Df := <0.00001617321837, 0.00001607066092>:
P__baro := <759.062, 759.062>:
PaCor := <94.82734101, 90.40928915>:
PvCor := <35.32630403, 35.55779803>:

MyProc:= proc(N, Hb, k, P50a, P50v, nu, Df, P__baro, PaCor, PvCor)
local p, P:= (2*p/(P50a+P50v))^N;
     Por*phi/4/(1-Por)/BP__length*(nu/Df)^(2/3)
     * int((1+1.34*Hb*N*P/k/p/(1+P)^2)^(2/3)/(P__baro - p), p= PvCor..PaCor)
end proc
:
MyProc~(N, Hb, k, P50a, P50v, nu, Df, P__baro, PaCor, PvCor);

Let E be the original set of edges. Then do

subsindets(E, symbol, String);

However, why not use GraphTheory:-RelabelVertices instead?

Making your corrected code into a procedure and getting matrix output can be done like this:

restart:
MyProc:= proc(pH, pO2, Tart, n)
local 
    P50p37:= 26.6*(10.^0.48)^(7.4-pH),
    S0:= (pO2/P50p37)^n
;
    [S0/(1+S0), (P50p37, pO2)*~(10.^0.024)^(Tart-37)]
end proc
:   
pH:= [7.398, 7.392]; pO2:= [121.6, 113.4]; Tart:= [32.5, 32.9];
n:= 2.7; 
Matrix(MyProc~(pH, pO2, Tart, n));
         [0.9836583446 20.78938475 94.82734141]
         [                                    ]
         [0.9799869417 21.39546041 90.40928912]

The multiplication that you want can be done by 

evalDG(5 * g5);

The Multiply command that you were trying to use comes from the LinearAlgebra package, and it has no connection to tensors. 

It took me less than 5 minutes to figure this out once I saw your worksheet. In the future, please be more forthcoming with any requested additional information.

You can't use square brackets [ ] or curly braces { } instead of parentheses (a.k.a. round brackets) for grouping and disamibiguation in ordinary algebraic expressions. Using [ ] or { } adds extra layers of meaning that you did not intend. It often takes several steps before an early use leads to an error, which is the case here. So, you need to start at the top and change all algebraic uses of [ ] and { }. 

You can handle products, quotients, powers, and incrementing the second argument like this:

evalindets[3](
    evalindets(
        evalindets(p, specfunc(T)^anything, `*`@op), 
        `*`, t-> ((A,B)-> add(A)*mul(B))(selectremove(type, [op](t), specfunc(T)))
    ),               
    specfunc(T), `+`, 2, applyop, 1 #This 1 is the increment.
);
  8 T(x, 8) + 8 T(x, 3) + 4 T(x, 6) + 11 T(x, 2) + 12 T(x, 4) + 7 T(x, 5)

The increment could be made -1 or anything else.

A multivariate truncated Taylor series can be obtained using mtaylor. Example (using your function):

f:= (x,y)-> x*y/(y-x*sqrt(y)-x^2):
mtaylor(f(x,y), [x,y]=~ [1,1], 6)

where [1,1] could be replaced by any point where f is analytic. Using instead

mtaylor(f(x,y), [x,y], 6)

implies that the desired expansion point is [0,0], which won't work for this particular function. 

It can be done by temporarily integerizing the exponents on epsilon. One such way is to substitute epsilon= epsilon^2 followed by combine(..., power). Immediately after defining A, do

A2:= collect(combine(subs(epsilon= epsilon^2, A), power), epsilon);
add(coeff(A2, epsilon, k)*epsilon^(k/2), k= 0..2);
 

I wouldn't use either way: The overload way leads to duplication of code, which causes big problems when code needs to be modified. The "traditional" way is bad because it's best to isolate the input processing from the main algorithm. In other words, within the main algorithm, you don't want to be taking branches due to trivial differences in the input type, because that distracts the reader from the finer mathematical points of the algorithm.

Maple provides a solution for this situation called "data type coercion" (help page ?coercion). If a procedure's argument is not the desired type, then another (trivial) procedure (that you write) will be called to convert it to the desired type, if possible. To use this, both the type name and conversion procedure name should begin with ~. Example:

restart
:
#Type-coercion procedure (name must begin ~): 
~set:= proc(A, T::type)
    if A::{list, thistype}(T) then 
        `if`(A::list, {A[]}, {A})
    else
        error "%1 not coercible to type set(%2)", A, T
    fi
end proc
:
foo:= proc(A::~set(`=`))
    A
end proc
:
foo({x=3, y=4}), foo(x=3), foo([x=3]);
                {x = 3, y = 4}, {x = 3}, {x = 3}
foo(3);
Error, (in foo) 3 not coercible to type set(=)

I think that @nm has the right idea for a "purely Maple" Answer. But I think that you should also learn the easy formula for it. Suppose that you have a decimal number of the form

x = n.d1d2...djr1r2...rk...

where n is nonnegative integer and d1, d2, ..., dj, r1, r2, ..., rk are digits (0--9) with the r1r2...rk being the repeating part. Then

x = n + (d1d2...dj + r1r2...rk/(10^k-1))/10^j

For example:

x:= 7.421232323;

r:= 7 + (421 + 23/99)/1000;
             r := 367351/49500

evalf(r);
             7.421232323

Here's a procedure for it:

F:= proc(A::Matrix, B::Matrix)
local A0;
    if upperbound(B)[2] < upperbound(A)[2] then return thisproc(B,A) fi;
    A0:= <A | Matrix(upperbound(A)[1], upperbound(B)[2] - upperbound(A)[2])>[.., ..-2];
    evalb({convert(A0, listlist)[]} = {convert(B[.., ..-2], listlist)[]})
end proc
:

Assuming that g is a binary variable (i.e.,  it can only be 0 or 1), the easiest way to deal with it is to solve separately for those cases, then take the maximum of the two solutions. It can be done like this:

params:= [
    Pp= 1600, delta= 8760, Cm= 104, x= 175200, Ep= 1.5, Rc= 8, v= 80, Ec= 0.5, Cs= 0.5,
    Csc= 0.6, Sc= 1500, Ch= 20, Er= 2, Rp= 100, Ec= 0.5, a= 100, b= 2.5, Do= 1500, Ig= 80,
    Ch1= 22, Ce= 6, Ces= 4, E(Pr)= 0.05, E(1/(1 - Pr))= 1.0536, E(Pr/(1 - Pr))= 0.0536
]:
TOTP := unapply(eval(ETPU(y,Sr), params), [Sr,y]):
TotP_g:= g-> Optimization:-Maximize(
    eval([TOTP(Sr,y), eval({300 <= D(g,Sr), D(g,Sr) <= 2000}, params)], :-g= g)[], 
    y= 200..400, Sr= 200..400
):
TotP_g(0);
    [94829.4555304706155, [Sr = 400., y = 317.129445700422]]

TotP_g(1);
  [114611.718817373156, [Sr = 400., y = 364.268415001416]]

Shouldn't the command point(C, a, -c) be point(C, a, -b)?

This is I think the most direct translation:

[ListTools:-Deal](
    [seq](
        `if`(x<0, x+2^16, x), 
         x= FileTools:-Binary:-ReadFile("Namefile", datatype= integer[2])
    ), 
    8
);

ListTools:-Deal combines Mathematica's Transpose@Partition into a single command. Deal was added in Maple 2021, so if you have older Maple, this'll need to be adjusted.

The output of FileTools:-Binary:-ReadFile is an Array. The [seq] converts that to a list. I only did that in the spirit of "direct" translation; an Array or Matrix would likely be better suited to your ultimate purpose. That would certainly be true if memory is an issue.

Maple's integer[2] format is signed 2-byte integers. The `if`(x < 0, x+2^16, x) converts them to unsigned.