In the following, assuming n::posint leads to the messy sum that you mentioned, and assuming n > 0 leads to a simple form in terms of Psi. But I can't convert those Psi to hypergeom.

restart:
evalindets(convert(1/(1+`x^n`), FormalPowerSeries), name, parse);
                      infinity           
                       -----             
                        \               k
                         )        k / n\ 
                        /     (-1)  \x / 
                       -----             
                       k = 0             

simplify(value(subsop(1= int(op(1, %), x= 0..1, continuous), %)))
    assuming n > 0;  #Assuming n::posint makes it too complicated.
                        /n + 1\      / 1 \
                     Psi|-----| - Psi|---|
                        \ 2 n /      \2 n/
                     ---------------------
                              2 n         

#Verify:
evalf(eval(%, n= 9) = Int(1/(1+x^9), x= 0..1));
                  0.9320304240 = 0.9320304242

If you use axesfont= [12, 12], the first 12 is assumed to be the name of a font, and the second 12 is the size of that font. You'd get exactly the same results using axesfont= [foo, 12]. A font specification in a plot command is always a list. Apparently, sometimes it's willing to ignore that the specified font 12 doesn't exist and just use some default font, and other times it's not willing to do that. Anyway, your usage of [12, 12] is simply nonsense.

I'm guessing that when you specified xk as the 2nd argument to solve, you were expecting it to infer that the solution variables were meant to be "all variables whose names begin xk." There is a (slightly advanced) way to do that, using a type called suffixed:

eq1 := xdot = (xk - xk1)/dt;
eq2 := xdot2 = (xk - 2*xk1 + xk2)/dt^2;
eq3 := c*xdot + k*xk + m*xdot2 = F;
eq:= {eq||(1..3)}; V:= indets(eq, suffixed(xk));
sol := solve(eq, V);

I assume that you mean that you want to find the index pairs [a,b] such that there exists an index x such that R[a,x] = 1 and S[x,b] = 1. Is that correct? In that case, you just need to multiply R.S and extract the index pairs of the nonzero entries of the product. This is amazingly easy in Maple. All that's needed is

[lhs]~(op(2, R.S))

Making an assignment to an essential predefined global variable such as constants is a little bit scary. I think that the preferred way to do what you want is to use the inert form of eval, named Eval:

mydiff:= Eval(diff(y(x,t), [x$3]), {x= a}) = V  # [...] is optional

That'll display the vertical evaluation bar in gray, which indicates the inertness. To display it in the usual blue:

InertForm:-Display(mydiff, inert= false)

I'd also prefer more automation, abstraction, and formalism for this. So, here it is:

restart:
Orth:= V-> [1,-1]*~V[[2,1]]:  #an orthogonal complement of a 2-vector
X:= [x,y](t):  P:= [r, theta](t):  #define variables
sys:= diff(X,t) =~ Orth(X);  #Cartesian-coordinate ODE system
#Use of a well-known and predefined transformation:
Tr:= X=~ changecoords(X, X, polar, P);
solve(eval(sys, Tr), diff({P[]}, t));

As you'd expect, PDEtools:-dchange works just as well for this, but I think it's overkill when working with a well-known predefined coordinate transformation.

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.