The eval command allows you to specify values for parameters that don't appear in the target. Because of this, there's no need to know the number of parameters that actually appear in the array as long as you know a superset of all possible parameters. Example:

PS:= [a, b, c]:  #parameter superset, as a list
Vals_1:= [1, 2, 3]:
Vals_2:= [4, 5, 6]:
Arr:= <1+a, b^2, a+b>:  #c doesn't appear in Arr
eval(Arr, PS=~ Vals_1);
eval(Arr, PS=~ Vals_2);
 

To facilitate the usage of step-by-step evaluation of inert expressions as a debugging tool, my program below implements it as follows:

1. Expressions can be copy-and-pasted directly from code. There's no need to use inert operators (%-operators) or to replace variables with values.
2. Values for the variables can be supplied in the call to the program.
3. Typesetting is used to highlight the subexpression being evaluated at each step, which'll appear in red.
4. Evaluation errors are trapped and displayed. The erroneous subexpression is retained in inert form and evaluation proceeds.
5. Subexpressions taking more than 1 second to evaluate are treated as in 4.
6. Evaluation proceeds in the natural order---innermost subexpressions first.

Here's the worksheet, which I can't display inline:

Download StepwiseEval.mw

And here's the code (this code cannot be entered as 2D-Input):

(*************************************************************************
StepwiseEval: 
    An appliable module for step-by-step evaluation of expressions,
    intended as a debugging tool

Input:
    X: 
        An expression entered as a string
    EV (optional): 
        Evaluations for variables in X. Any equation (or set or list of
        equations) allowed as the 2nd argument to eval is allowed.

Displayed output:
    1. A typeset version of each step of evaluation with the part that'll
    be next evaluated highlighted in red. (Multiple occurences of 
    identical subexpressions are processed simultaneously.)
    2. The final value.
    3. Errors during evaluation are ignored and the error messages are
    printf'd. Each subexpression evaluation is timelimit'd to 1 second.

Returned output:
    A Record with 3 fields:
        'final': The final value
        'raw_forms': A vector of each step of evaluation as an ordinary
            Maple expression
        'displays': A vector of the typeset and highlighted versions of
            the raw_forms and string forms of any error messages

**************************************************************************)

StepwiseEval:= module()
option `Author: Carl Love <carl.j.love@gmail.com> 2021-Nov-11`;
export
    #Extract underlying "name", i.e. symbol, of a function. The loop is
    #needed to handle cases such as %f[a](b)(c).
    RootSymb:= proc(f::function)::symbol;
    option cache;
    local r:= f; 
        do r:= op(0,r) until r::symbol 
    end proc,

    #type for functions whose names begin `%`:
    ModuleLoad:= proc()
        if not TypeTools:-Exists(%Func) then
            TypeTools:-AddType(
                %Func,
                proc(f) 
                option cache;
                    #..1 instead of 1 to handle null function ``.
                    f::function and (""||(RootSymb(f)))[..1] = "%"
                end proc
            )
        fi;
        return
    end proc, 

    ModuleApply:= proc(X::string, EV::{thistype, list, set}(`=`):= {})
    local
        Value:= proc(f::%Func) 
        option remember;
        local s:= RootSymb(f);
            try 
                timelimit(1, subs[eval](s= substring(s, 2..), f)) 
            catch:
                printf(
                    "%s",
                    (dres,= sprintf(
                        "Error bypassed while evaluating %a: %s\n", 
                         f,
                         StringTools:-FormatMessage(lastexception)
                    ))[-1]
                );
                f  #Return original inert function.
            end try
        end proc,
  
        #This subsindets is needed because Parse treats all `-` as unary;
        #we need binary %- for this program. Specifically, all 
        #occurences of a %+ (-1)*b are converted to a %- b, and the
        #equivalent is done for %+ with 3 or more arguemnts.
        old:= subsindets(
            eval(InertForm:-Parse(X), EV), 
            specfunc(`%+`), 
            curry(
                foldl, 
                (a,b)-> 
                    if b::negative or b::`*` and op(1,b) = -1 then 
                        a %- (-b) 
                    else 
                        a %+ b
                    fi
            )@op
        ), 
        res:= <old>, dres:= rtable(1..0, subtype= Vector[column]),
        temp, f, df,
        grey:= "#909090", red:= "#FF0000" #color codes used by Typesetting
    ;
        do
            #Sorting by length causes inner functions to be done 1st.
            for f in sort([indets((temp:= old), %Func)[]], length) do
                temp:= eval(temp, f= Value(f))
            until temp <> old;
            if temp = old then
                print(old); 
                return Record(
                    "final"= old, "raw_forms"= res[..-2], "displays"= dres
                ) 
            fi;
            print(
                (dres,= eval(
                    InertForm:-Display(old), 
                    (df:= InertForm:-Display(f))= subs(grey= red, df)
                ))[-1]
            );
            res,= (old:= temp)
        od
        end proc
;
    ModuleLoad()
end module
:    

#Usage example:
#This is the same expression that Acer used above except that
#   - it's a string,
#   - no % operators are used,
#   - some numbers have been replaced by variables.

S:= "A/(1 - abs(B^2 - (C - 1 - 1^2)^2)) + abs(E)/abs((x - x) - 2^2)^2"
:
#The 2nd argument below is variable values to make the expression identical to Acer's
#example:
Res:= StepwiseEval(S, {A= 12, B= 2, C= 3, E= -7}):

The typeset and prettyprinted output can't be displayed on MaplePrimes. You can see it in the attached worksheet.

To do what you're trying to do above, your for loop can be replaced by a single elementwise operation. All that's needed is

aVectorList:= [<1,2>, <3,2>];
results:= rtable~(1..4, aVectorList)

The 2nd line could also be

results:= Vector~(4, aVectorList)

although Vector is less efficient than rtable (since Vector is a high-level procedure with extensive argument processing whereas rtable is built-in).

Using rtable with an initializing vector (or rtable~ with a container of initializing vectors) as above is semantically equivalent to (although more efficient than) using copy or LinearAlgebra:-Copy. Thus, the mutability issue is avoided. 

The trailing 0s in the vectors come from the fill option to either rtable or Vector, which defaults to fill= 0. Thus this method avoids the potential issue of stale entries not being zeroed, as mentioned by acer.

Elementwise operation is often the most-convenient way to avoid the extreme inefficiency of building a list by iteratively adding elements to an existing list.

Your ODE can be easily solved symbolically by the methods of first-year calculus: just simple integration. The solution is an elementary function. All those exp functions are just obfuscation because they don't depend on y. All that matters is

ode:= diff(theta(y), y$2) = -(A*sinh(M*y)+B*cosh(M*y))^2;
int(rhs(ode), y);
S:= int(%+C1, y)+C2;
Cs:= solve({eval(S, y=-sigma)=0, eval(S, y=sigma)=1}, {C1,C2});
combine(simplify(eval(S, Cs)));

Then substitute the appropriate expressions for A and B.

If you use dsolve(..., numeric) for this (or any) BVP, it will use a finite difference method automatically.

Another way to avoid "deep" evaluation is to do it in a procedure:

a:= 1:
save a, "A.m":
restart:
NewNames:= proc(com::string)
local OldNames:= {anames('user')};
    parse(com, statement);
    {anames('user')} minus OldNames
end proc
:
b:= 2:
NewNames("read `A.m`;");
             {a}

You have unbalanced parentheses in ODE11, ODE12, and ODE13.

I see this as a bug (or flaw) in the design of DataFrame, but given that that design has been implemented and that that implementation has been released for several years, it wouldn't be easy to fix it and maintain backward compatability.

@vv wrote: Probably such numerical labels should not be accepted.

I disagree. Here is how I would've designed it: Allow anything to be labels, but if labels are used then those labels rather than row and column numbers should be the only way to index. While it would be difficult to change DataFrame itself in this manner (because of the backward compatibility), it would be very easy to write a wrapper for it that implemented this idea. Since a DataFrame is itself just a wrapper for a matrix, it'd still be trivial for a programmer to index the underlying structure by rows and columns.

I worked on large databases with significant user interfaces (point-of-sale systems, inventory, etc.) for several years in the late '80s and early '90s. I learned the following the hard way: There is no good reason that the end user should ever see, let alone refer to, the actual record numbers (or any other field that is used as the primary key). Indeed, if you allow them to see those keys, several problems will gradually arise:

1. Gradually (and perhaps at first subconsciously) they begin to ascribe significance or meaning to certain numbers or keys or ranges thereof. That's just the way the human mind works.
2. Once that happens, occasions will arise where they will want to change some of those keys. Real-life example: A customer's status changes from retail to contractor. The user wants to change the customer's number in whatever way indicates to them that that's a contractor, even though I've explained to them that it'd be much less work for me (and cost for them) to just add a 1-bit field for contractor status.
3. Doing that can be a major project because those numbers (or primary keys) need to be linked to many other files in the overall database.

The way around this is to use a secondary key, which will seem to the end user to be the primary key, yet they can do however they please with it. But the secondary key only exists in that one file (and its index file). All linkages to other files are via the primary key. In the case of DataFrame, the labels are those secondary keys.

One way:

ST:= table((parse@lhs=rhs)~(S)):
MS:= [mu, sigma]:
MS =~ index~(ST[2], MS); rhs~(%)[];

In cases where the known function f(x) has certain relatively simple forms (such as polynomial), a usable solution for unknown g(x) can be obtained like this:

restart:
f:= x-> 2*x+5:
h:= f*g:
h2:= (D@@2)(h)(x);
                 h2 := 4 D(g)(x) + (2 x + 5) @@(D, 2)(g)(x)

#We need to "hide" the symbolic derivatives because they 
#will confuse the solve command:
S:= {D(g)(x)= g1, (D@@2)(g)(x)= g2}: 
h2i:= subs(S, h2);
                         h2i := 4 g1 + (2 x + 5) g2

sol:= solve({h2i < 0, g1 > 0, g2 < 0}, x, parametric);

          sol :=  / [[  4 g1 + 5 g2    ]]                         
                  | [[- ----------- < x]]      And(g2 < 0, 0 < g1)
                 <  [[     2 g2        ]]                         
                  |                                               
                  \          []                     otherwise     

#Thus, the x-interval(s) of concavity for h are where this inequality
#is true:
Concavity_condition:= op(2, sol)[][];
                                           4 g1 + 5 g2    
                  Concavity_condition := - ----------- < x
                                              2 g2        

#Since we know g2 < 0, that can be simplified to
Concavity_condition2:= %f(x)*g2 + 4*g1 < 0
                 Concavity_condition2 := %f(x) g2 + 4 g1 < 0

#In this case, I did that simplification by hand. 
#Note that I used % to make f(x) inert.

If you can make that table into a 3-column Matrix (I'll call it A), then do

plots:-pointplot3d(A, connect)

Or if you can make three Vectors X, Y, Z, then do

plots:-pointplot3d(X, Y, Z, connect)

Or if you can make the table into a list L of 3-element sublists, then do

plots:-pointplot3d(L, connect)

Like this:

restart:
interface(prettyprint= 1):

#Your original input:
ode:= diff(y(x),x) = 2*(2*y(x)-x)/(x+y(x)):
ic:= y(0)=2:

#Transform to inverse function:
ode1:= 1/diff(x(y),y) = subs(y(x)= y, x= x(y), rhs(ode));
ic1:= (op@lhs=x@rhs)(ic);

#dsolve(implicit), then transform back:
sol:= subs(x(y)=x, y= y(x), dsolve({ode1,ic1}, 'implicit'));

                                 1       2 (2 y - x(y))
                      ode1 := -------- = --------------
                               d            x(y) + y   
                              --- x(y)                 
                               dy
                      
                               ic1 := 0 = x(2)

               /2 x - y(x)\       /x - y(x)\                              
   sol := -3 ln|----------| + 2 ln|--------| - ln(y(x)) + I Pi + ln(2) = 0
               \   y(x)   /       \  y(x)  /                              

#A little bit of algebra shows that that's equivalent to your
#"book solution":
exp(combine(lhs(sol), symbolic)) = exp(rhs(sol));
                                           2    
                               2 (x - y(x))     
                             - ------------- = 1
                                           3    
                               (2 x - y(x))     

Here's a procedure for it. It counts all the indices in one pass:

Count:= (C::{name, name^posint, `*`({name, name^posint})})->
    (Statistics:-Tally@op~@[op]@indets)(
        subsindets(
            C, And(name, indexed)^posint,
            e-> `*`(
                'convert(op([1,0], e), `local`)[(op@op)(1,e)]' 
                $ op(2,e)
            )
        ),
        indexed
    )
:
C:=x*u[]*u[1,2,2]^3*u[1,1,2]:
Count(C);
                         [1 = 5, 2 = 7]

As nm and acer have explained, the difference that you're seeing is not due to a difference between map and ~ but due to the difference between matrix/vector multiplication in the old linalg package vs the modern A.B. However, there are two differences between map and ~ that average users should be aware of:

1. ~ can only map over lists, sets, rtables (which includes Vectors, Matrixes, and Arrays), and tables (I'll call these the neutral containers for the rest of this post); whereas map can map over almost any composite structure (functions, polynomials, etc.). Function arguments that aren't one of those types get treated as singletons. Example 1a: Let L be a list of pairs with each pair being a list; then op~(2, L) returns the same as map2(op, 2, L). Example 1b: The list of prime factors of any integer n is returned by op~(1, ifactors(n)[2]). In both these examples, the first argument of op, 1 or 2, is being ignored by ~ (and still seen by op).
2. map only maps over one object; whereas ~ maps simultaneously over all neutral containers that appear as arguments as long as their sizes (and sometimes shapes) conforms. Example 2: f~([1,2], 3, {a, b}).

Any plots of the same number of dimensions (either 2 or 3) can be combined into a single plot with the command plots:-display. The plots don't need to have anything in common other than the number of dimensions. So, you could do this:

p1:= plots:-arrow(...);
p2:= plots:-arrow(...);
p3:= plots:-arrow(...);
plots:-display(p1, p2, p3, ...);

where the ... in the arrow commands represents exactly the arguments in your posted worksheet, and the ... in the display command represents possible other options that would apply to the plot as a whole (such as scaling= constrained). The above could also be done as

plots:-display(
    plots:-arrow(...),
    plots:-arrow(...),
    plots:-arrow(...),
    ...
);

I didn't see any vector field defined in your worksheet; however, if there had been, you could use one the several commands for plotting vector fields (such as plots:-fieldplot3d) and include that in the display command.

All the information is in the paper, but the portion that you showed only gives this obviously periodic function:

x:= t-> exp((cos(t)-cos(T))/4);

They say T=Pi/2, t>T, and gamma=0.4 (although your screenshot doesn't show where gamma is used). We can plot x like this:

T:= Pi/2:
plot(x(t), t= 1..24, view= [1..24, -0.2..1.8]);