@sursumCorda 

I was writting an explanation but was in the end beatten by @acer

Without any excursion into the third dimension:

pt := [
        seq([u, 5*(u-1)], u in [seq](1..2, 0.01))
        , seq([u, 5], u in [seq](2..1, -0.01))
        , seq([1, v], v in [seq](5..0, -0.01))
      ]:



S  := display(polygon(pt, color=black)):
plots:-display(
  transform((s,t)->[s^2*sqrt(t)*cos(t), s^2*sin(t)])(S)
  , size=[500,"golden"]
  ,axes=boxed
);

MyExplanation.mw

Here is a way

Download Numeric_mmcdara.mw

A step by step demonstration of the content of the math book

If your goal is to implement in Maple what is written in your math book as closely as possible, here is a solution.
 

Download StepByStep_Procedure.mw


If you are in search of different methods to get the solution, here is a first one

The code below founds 4 solutions from which only one corresponds to a strictly positive volume.
(replace L by 6)
 

restart:
v := a*b*c: 
p := a+b+c:
# with Lagrange multiplier lambda 
U := {a, b, c, lambda}: 
diff~(v-lambda*(p-L), U): 
sols := solve(%, U); 
         {a = L, b = 0, c = 0, lambda = 0}, 

           {a = 0, b = 0, c = L, lambda = 0}, 

           {a = 0, b = L, c = 0, lambda = 0}, 

            /    1        1        1             1  2\ 
           { a = - L, b = - L, c = - L, lambda = - L  }
            \    3        3        3             9   / 

sol  := remove(has, [sols], lambda=0)[]
           /    1        1        1             1  2\ 
          { a = - L, b = - L, c = - L, lambda = - L  }
           \    3        3        3             9   / 

and a second one based on a derivative-free demonstration: 

# Volume of the cube a=b=c=L/3

Vc := (L/3)^3:

# Volume of the parallelepiped 
# a = L/3-delta__a
# b = L/3-delta__b
# c = L/3-delta__c


Vp := (L/3-delta__a)*(L/3-delta__b)*(L/3-delta__c):

# Variation of volume regarding to Vc

delta__v := Vp - Vc:

# Let's account for the condition delta__a + delta__b + delta__c = 0

delta__v := expand( eval(delta__v, delta__c=-delta__a-delta__b) );

# Let's arrange the expression

map(factor, collect(delta__v, L));

         1 /        2                               2\  
       - - \delta__a  + delta__a delta__b + delta__b / L
         3                                              

                    2                             2
          + delta__a  delta__b + delta__a delta__b 

# The second order perturbation is obviously negative
# for delta__a^2+delta__a*delta__b+delta__b^2 > 0
#
# Then the variation of the volume of the parallelipiped at 
# point a=b=c=L/3 is strictly negative unless delta__a, delta__b=0
# and thus delta__c are both null.
#
# The cube is the parallelipiped with the highest volume given 
# the value of L.

I'm not sure I really understood your question.
So if my answer is off the mark, forget it quickly:

(the flag see enables removing the outputs)

Download SomethingLikeThat.mw

Is it this that you want to achieve?

Eq := I__C*H(t)-m*rho*rho*H(t)
lprint(Eq);
Physics:-`*`(I__C, H(t))-m*Physics:-`*`(Physics:-`^`(rho, 2), H(t))

eval(Eq, `*`=:-`.`):
lprint(%);
I__C . H(t)-m*(Physics:-`^`(rho, 2) . H(t))

I do not understand why you talk about the Physics package.

Is it this package that you really want to use (seems to me a little be complicated for the question you ask).
If you want something else why not

restart
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 

   21 2015 Build ID 1097895
with(VectorCalculus):
v := <c, d>;  # the display is smarter in a worksheet
           v := c*e[x]+d*e[y]
position := <p, q > + int(v, s)
           position := (c*s+p)*e[x]+(d*s+q)*e[y]

Without any package;

restart
v := <c, d>;
p := t -> <x(t), y(t)>:
         
position := eval(p(t), dsolve({ diff~(p(t), t)=v, p(0)=<alpha, beta> }))
          

Concerning "_i"

You can change _i by _WhatYouWant and get exactly the same kind of results (keep in mind what I told you yb a prevo-ious question about names beginning with an underscore https://www.mapleprimes.com/questions/235461-What-Is-The-Type--obtained-For)
_i has absolutely nothing particular in general:

about(Pi)
Pi:
  is assumed to be: Pi
  a known property having {irrational} as immediate parents
  and {BottomProp} as immediate children.
  mutually exclusive with {0, 1, integer, prime, composite, fraction, rational}

about(_i)
_i:
  nothing known about this object

When you use the Physics package things may change.
See  < Overview of the Physics[Vectors] Subpackage > :

Regarding projected vectors, the Vectors subpackage is designed to work only with cartesian, cylindrical and spherical orthonormal basis and the related systems of coordinates (see Identify), according to
           (_i, _j, _k)
          ...
and further
NOTE: these variables x, y, z, rho, phi, r and theta, as well as _i, _j, _k, _rho, _phi, _r and _theta, respectively used to represent the coordinates and the unit vectors, are automatically protected when the Physics[Vectors] subpackage is loaded.


Concerning "D"

# here f is explicitely defined as a function a a single argument u
f := u -> u*v:
D(f)
                          u -> v

# nothing but the product rule you learnt at school
# https://en.wikipedia.org/wiki/Product_rule
g := u*v:
D(g)
                        D(u) v + u D(v)

Finally, the correct syntax is not D(f(u)) but

D(f)(u)
                               v
D(f)(t)
                               v

Look to D help page to see how to use D correctly
When tou write D(f(u)) the first step is the evaluation of f(u), so it is  D(u*v) that you compute, and thus the result you got.

Question 1
Given two vectors A and B of same type (both row vectors or both column vectors) and same "length", and given two numbers p and q then

C := p*A + q*B

buids a vector C of the same type and length whose element C[i] = p*A[i]+q*B[i]

A := <1 | 2>:
B := <2 | 1>:
C := 2*A + 3*B;
           C :=  [ 8 7 ]

For other componentwise operations on vectors or matrices use ~:

E := A *~ B
            E := Vector[row](2, {(1) = 2, (2) = 2})
F := cos~(A) *~ exp~(B)
            F := Vector[row](2, {(1) = cos(1)*exp(2), (2) = cos(2)*exp(1)})

Question 2
About _names the things are very clear: never ever use such names.
If you are fan with underscores as prefix, just double them:

__a := 3:
__b := 4:
__c := __a + __b;
                               7

This will spare you a lot of problems.
Even if a single underscore can lead to no error at all, your are never sure of that:

# all is fine here
_toto := 3:
_toto^2
                               9

Question 3
Look to the type help page and click on +:

• type/`+` - check for an expression of type `+`
• type/`*` - check for an expression of type `*`
• type/`^` - check for an expression of type `^`

The fact c:=a+b is of  `+` type comes from the way an algebraic operation (here a sum) is represented by a syntactic graph:

c := a+b:
whattype(c)
                               +
ToInert(c)
         _Inert_SUM(_Inert_NAME("a"), _Inert_NAME("b"))

... unless you have a problem with your version.
(BTW your worksheet doesn't contain any plot command: wher do your plots come from ???)

In the RootVector help page search for the string plot and do what it's said:

Download PlotRootedVector_mmcdara.mw

I guess the ambuity comes from the use of the names "create" and "indexed".

Let's take the case of a procedure through two examples. 
In the first one procedure f assigns to a variable a the value 1, but Maple returns a warning telling you it considers a as a local variable. In the second example a is explicitely declared local

f := proc()
a := 1:
end proc
Warning, `a` is implicitly declared local to procedure `f`


interface(warnlevel=0): 
a:=2: 
f := proc() 
a := 1: 
end proc: 
f(); 
a;
                               1
                               2

f := proc()
  local a:
  a := 1:
end proc

a := 'a'; # 
f(); 
# a is still unknown out of procedure f
a; 
# a way to assign 
assign(a=f()): 
a;
                               a
                               1
                               a
# and another one
print():
a:='a'; 
a:=f(): 
a;              
                               a
                               1
                               a

From help page < local > ... If present, specifies the local variables reserved for use by the procedure...
From help page < assign >... assign(a, B) and assign(a = B) commands make the assignment a := B and return NULL.

assign associates a name (a) with the value of an expression (B).

About the lines a:='a' the mechanism of this clearing operation is illustrated here

restart:
anames(user);
a := 1;
anames(user);
                               1
                               a
a := 'a':
anames(user);  # nothing returned, a no longer exists

The use of local (or global too)  enable declaring names, not really (IMO) creating variables ; a name can be associated with an expression and it's common to say this name is a variable
(
I don't know is this reference can be found freely on the web "Introduction to Maple, André Heck, Springer-Verlag, 1993" ?
If it is so, from page 63: "With the procedure assign you can get the status of individual variables".
Then assign, or :=, both give to a name the status of variable ; I guess one can say that assign creates variables.
)

According to the book "Maple V, language Reference Manual, Char et al., 1991", the term variable is used quite lately in it, mainly in refereces to the previous clearing process and he use of the global statement.

global acts basically as local does, excepted that the use of global a; makes a to be known at all levels in your worksheet:

restart:

f := proc()
global a:
a:=1:
end proc:

anames(user);
f():
a;
                               f
                               1

From verion 2019 (I think) it is no monger possible to use the statement global at the upper level of the worksheet, that is to write something like this

restart:
global a:
                               f
                               1

So let me rewrite your question as  "Is it at all possible to assign indexed names"
The answer is yes. From < indexed > help page:

• An indexed expression of the form name[expression sequence] represents a selection operation.
• An indexed name is a valid name.  In particular, since A[1,2,3] is a valid name, so are A[1,2,3][x,y] and A[1,2,3][x,y][2,1].
• The use of the indexed name b[1] does not imply that b is an array.  For instance, if b is an unassigned name, a := b[1] + b[2] + b[1000] simply forms a sum of three indexed names.
• If an assignment is made to name, then name[expression sequence] represents a selection from name using the index sequence expression sequence.  For more information, see selection.

Here are a few valid names:

• R1, R2, R3, ...
• R[1], R[2], R[3], ...

Only the last ones are indexed names:

restart:
R[1] := 1:
R1   := 2:
R[1], R1;  # R[1] is not R1 !!!
print("--------------------------"):

restart:
type(R[1], indexed);
type(R1, indexed);
                              1, 2
                  "--------------------------"
                              true
                             false

An important thing (see bove) is that "The use of the indexed name b[1] does not imply that b is an array" (nor a vector).
The index is not necessarilly an integer:

restart:
type(R[p], indexed);
# but...
type(R[p], indexed[integer]);
                              true
                             false

Be carefull, some names which appear as indexed names are nor indexed names:

restart:
L := [R[1], R1, R__1];
lprint(L);          # look to < lprint > help page
type~(L, indexed);  # look to < ~ > (operators, elementwise (~) ) help page
                        [R[1], R1, R__1]
[R[1], R1, R__1]
                      [true, false, false]

How to assign indexed and non indexed names?
What won't work:

seq(R[n], n=1..3);
f := proc()
  local seq(R[n], n=1..3)
end proc:
                        R[1], R[2], R[3]
Error, `(` unexpected

f := proc()
  local R[1], R[2], R[3]:
end proc;

Error, `[` unexpected

seq(R||n, n=1..3);
f := proc()
  local seq(R||n, n=1..3)
end proc:
                           R1, R2, R3
Error, `(` unexpected

seq(R__||n, n=1..3);
f := proc()
  local seq(R__||n, n=1..3)
end proc:
                        R__1, R__2, R__3
Error, `(` unexpected

And what works

f := proc()
  local R1, R2, R3:
end proc;

f := proc()
  local R__1, R__2, R__3:
end proc;

but neither R1 nor R__1 are indexed names.

If you know the number n of indexed names you want to create 

restart
f := proc(n)
  local R := Vector(n, symbol=R):
  local i:
  for i from 1 to n do R[i]:=i end do:   return R; end proc:

res := f(3); 
            Vector(3, {(1) = R[1], (2) = R[2], (3) = R[3]})

assign(seq(R[i]=res[i], i=1..3)):  # multiple assignment

# check for the values of R[1], R[2] and R[3]
'R[1]' = R[1];
'R[2]' = R[2];
'R[3]' = R[3];
                            R[1] = 1
                            R[2] = 2
                            R[3] = 3

Finally, and it is how I understant you question, you can associate indexed names with expressions at the upper level while doing this

restart:
assign(seq(R[i]=t^(i-1)*cos((i-1)*Pi/2), i=1..3)):

'R[1]' = R[1];
'R[2]' = R[2];
'R[3]' = R[3];
                            R[1] = 1
                            R[2] = 0
                                    2
                           R[3] = -t 

and, more specifically:

restart:

MyFiles := FileTools:-ListDirectory(cat(currentdir()))[1..3];
        MyFiles := ["EgyptianFraction.mw", "matrix_inverse_mmcdara.mw", "Formal_ODE_Solution.mw"]

assign(seq(R[i]=MyFiles[i], i=1..3)):

'R[1]' = R[1];
'R[2]' = R[2];
'R[3]' = R[3];

                  R[1] = "EgyptianFraction.mw"
               R[2] = "matrix_inverse_mmcdara.mw"
                R[3] = "Formal_ODE_Solution.mw"

Since I am not familiar enough with the Physics package, here is an workaround, for what it is worth.

e := expand(lhs(Newton_generalized_MatForm)):
select(has, [op(convert(e, D))], :-D(`ϰ_`)(t)):
convert(coeff(eval(add(%), :-D(`ϰ_`)(t)=J), J), diff)

Nevertheless I can't advise you enough to await a better solution ( @ecterrab ? )

If you want to store several Maple-generated matrices converted into Python code, I propose you this:

restart:
# just a notional example 
with(LinearAlgebra):
t := table([seq(i=RandomMatrix(4, 4), i=1..100)]):

# redirect the default output (terminal) to a file named "python.txt"
f := cat(currentdir(),"/Desktop/python.txt"): 
writeto(f):
for i in sort([indices(t, nolist)]) do
  CodeGeneration:-Python(t[i])
end do:

# reset the default output for subsequent work, if any
writeto(terminal);

python.txt

If you want to read from Maple the Python matrices the file above contains, do this

restart
f    := cat(currentdir(),"/Desktop/python.txt"):
T    := table([]):
line := readline(f):
i    := 1:
while line <> 0 do
   if line <> "> writeto(terminal);" then
     T[i] := convert(parse(StringTools:-Split(StringTools:-Split(line, "(")[-1], ")")[1]), Matrix):
   end if:
   line := readline(f):
   i    := i+1:
end do:
eval(T)

PythonMatrices_mmcdara.mw

It is always possible to rename the vertices as you want:

Download ExtensionField_mmcdara.mw

Generating automatically the vertex names is more complex while not impossible.
Note that non-breaking spaces can also be added in the vertex names for a better rendering and that their colors can be adjusted the same way.

Executing showstat(GroupTheory:-DrawSubgroupLattice) reveals DrawSubgroupLattice offers more capabilities then the corresponding help page (Maple 2015.2) describes (look to SubgroupLattice for more help).

  {
    assignlattice::name := undefined, 
    center := undefined, 
    derived := false, 
    directed::truefalse := false, 
    [edgecolor, edgecolour] := "grey", 
    highlight::{object, set, list({object, set, {object, set, list({object, set})} = anything}), {object, set, list({object, set})} = anything} := [], 
    indices::truefalse := false, 
    labels::identical(:-integers,:-none,:-zuppos,:-ids) := :-integers, 
    normal := true, 
    output::{list(identical(:-graph,:-plot)), 
    identical(:-graph,:-plot)} := :-plot, 
    title := undefined, 
    titlefont := undefined, 
    [vertexcolor, vertexcolour] := "white"
  }

In particular, if you want more flexibility in customizing the graph:

with(GraphTheory):
Gr := DrawSubgroupLattice(g,  'output'=':-graph');
# Then you can use all the features GraphTheory proposes,
# for instance:
Edges(Gr);
GraphTheory:-Vertices(Gr);

Other way:

L  := SubgroupLattice(g);
Gr := convert(L, 'graph')

Even if I guess you have already thought about it

map(D[1], z(u,v)):
type(%, Vector);
      true

diff_mmcdara.mw

For fun (maybe...)

Your polyhedron is a half truncated octahedron 
https://en.wikipedia.org/wiki/Truncated_octahedron
https://mathcurve.com/polyedres/octaedre_tronque/octaedre_tronque.shtml  (in french)

The title refers to the equation its boundary verifies (last reference above)

HalfTruncatedOctahedron.mw

This worksheet also contains a few lines to display all the poyhedra that plots:-polyhedraplot knows.
Could be useful