The entire thing can be done without constructing any arrays (or any other containers of graphs), without doing any counting, and without doing any indexing. Indeed, it can even be done without constructing any graphs; it's better to work directly from the adjacency matrices.

A graph's edge set is in both Es and Fs  iff
vertex n doesn't appear in the edge set  iff
the degree of vertex n is 0  iff
the nth row of the graph's adjacency matrix is all 0.
Thus, the following code does what you want: 

n:= 4:
currentdir("C:/Users/carlj/desktop"): #Change this line!!
L:= GraphTheory:-NonIsomorphicGraphs(
    n, output= iterator, outputform= adjacency, selectform= adjacency,
    select= (A-> evalb(ArrayTools:-AnyNonZeros(A[-1]) = 1))
):
Text:= 
    "\ntext1 text2 text3\ntext2 text2 text3\ntext3 text3 text3\n"
    "text4 text4 text4\ntext5 text5 text5\n\n"
:
fp:= FileTools:-Text:-Open("Graphs.txt", create, overwrite):
while (A:= L()) <> FAIL do 
    fprintf(fp, "G=%a"||Text, {lhs}~(op(2,A))) 
od:
FileTools:-Text:-Close(fp): 

You must change that currentdir command to something relevant to your system.

I'm on my phone right now, so I can't test this Answer in Maple. I will test when I get home.

The result that you got is not numeric, and it has no relation to Runge-Kutta or any other numeric method.

Answer 1: The solution returned by dsolve contains unevaluated integrals. This is often the dsolve response to a non-numeric BVP (initial conditions specified at two x-values) even when the underlying ODE is trivial, as this one is. These integrals themselves appear to be trivial. By applying the value command to dsolve's result, you should be able to force it to do the integrals. 

The fact that the integral signs are displayed in gray rather than blue indicates that these integrals have a possibility of further evaluation using value. That further evaluation doesn't always produce useful results, but in this case the integrals are trivial so it does.

Answer 2: Remove the second initial condition, Uy(L)=0, from the dsolve command. The returned result should have one constant of integration (with a name like _C1). Solve for that constant by applying the second initial condition. 

You could put the call to the inner procedure in an if statement, like this:

foo:= proc(n::integer)
    if
        proc()
            if n<0 then return FAIL fi;
            print("doing something else inner")
        end proc()
        <> ()
    then return 
    fi;
    print("doing something else outer")
end proc
:

In this case, the FAIL could be replaced by anything other than NULL (or its equivalents). Note that all executed procedures either end with an error or have a return value, possibly NULL, regardless of whether they've executed an error or a return statement. Thus, the invocation of any procedure can be placed in the statement

if proc() ... end proc <> () then ... fi

So there's not much point in invoking an anonymous procedure without incorporating its return value in an expression.

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]