In your Question, you present the following vector of equations (or expressions implicitly equated to 0) to be solved for the w variables. The equations are linear in those variables.

equ:= < 
    0.85*((phi__cs*beta__s*f__con) . (D__pile*w__1)) + <-5/6*H - 5/8*P>, 
    0.85*((phi__cs*beta__s*f__con) . (D__pile*w__2)) + <5/6*H - 5/8*P>, 
    0.85*((phi__cs*beta__s*f__con) . (D__pile*w__3)) + <1/2*H - 3/8*P>
>;

The form of these expressions is quite awkward. Before proceeding, it's mandatory to correct these 2 things:
1. `.` multiplication must be changed to `*`.
2. Internal vectors must be converted to scalars.

There's another optional correction that I'd like to make, because it makes the output nicer:
3. Convert decimal numbers to exact fractions.

All three corrections can be done automatically, regardless of the number of entries:

eqs:= [evalindets](
    evalindets(
        convert(equ, rational), #Convert 0.85 to 17/20.
        specop(`.`), `*`@op #Correct improper multiplication.
    ),
    rtable, seq #Convert internal vectors to scalars; convert outer vector/matrix to list
);

Your equations are of course trivial to solve individually. However, I suspect that in general you will want to have equations that contain more than one of the w variables each. The code below will also handle that more general case (as long as they're still linear in those variables).

LA:= LinearAlgebra:
W:= indets(eqs, suffixed(w__))[]; #Extract decision variables.
<W> =~ LA:-LinearSolve(LA:-GenerateMatrix(eqs, [W]));

You say that the awkward format of your original equ is due to it being the output of some other program (or equation). But there is no good reason why that program should produce its output in that form. You should change that other program or get its author to change it. 

Constructing matrices from "blocks" (or submatrices) and permuting the rows and columns of matrices is much easier than anything shown in this thread so far, and it can be done using no packages (no LinearAlgebra, no ListTools, no linalg), no loops, and no index variables (such as the i, j, and k that you used).

B1:= 24*<9, -12; -12, 16>; #2x2 block
Z:= Matrix(2,2); #2x2 zero matrix
#3x3 arrangement of 2x2 blocks to produce 6x6 matrix:
K1:= <B1, -B1, Z; -B1, B1, Z; Z, Z, Z>;
K2:= K1[[5.., ..4]$2]; #Apply permutation [5,6,1,2,3,4] to rows and columns.
B3:= <500, 0; 0, 0>;
K3:= <B3, Z, -B3; Z, Z, Z; -B3, Z, B3>;

#This creates a flattened form of your KG:
KG__flat:= <K||(1..3)>;

#If for some strange reason you actually want the non-flat
#KG that you originally had:
KG:= rtable(1..3, 1..1, [K||(1..3)], subtype= Matrix);

Here are some answers to some of your 3 additional questions:

1. How to flatten matrices?
Rhetorical answer: It's much easier to create a flat matrix in the first place (for example, as shown above) than it is to flatten a matrix whose entries are themselves matrices or vectors. So, is there some reason why you can't create them flat in the first place? I'm not saying that flattening is impossible, just that it'd be better to avoid it.

2. How to count the rows or columns of a matrix?
Answer: There are many ways to do it. Here are three ways to count the rows of a matrix M:

[op](1, M)[1];  [upperbound](M)[1];  rhs([rtable_dims](M)[1]);

To count the columns, replace [1] by [2] in any of the above.

3. How to delete rows and columns from a displayed matrix to facilitate the entering of new entries in 2D Input?
Non-answer: I don't know how to do things in 2D Input. However, it's trivial to delete rows or columns, or make them blank, or just re-enter them in 1D input. 

See ?ProgrammingGuide,Appendix1. In particular, go to "Internal Representation" => "Internal Representations of Data Types" => "BINARY".

BINARY is one of the 62 fundamental data types in Maple (63 if you count "temporary"). Many of those data types cannot have a pointer pointing at them that can exist at the user level. In other words, they cannot appear as the right side of an assignment statement in Maple code, nor can they be otherwise manipulated in Maple expressions; they can only be processed in internal kernel code. BINARY is one of those types.

kernelopts(dagtag = BINARY);
             56

pointto(assemble(56));
Error, object at address is binary

The call to p1 inside p2 should be 

p1(a, b, ':-debug'= debug)

The ':-debug' refers to the global keyword; the debug refers to p2's parameter. Putting a parameter in quotes, as your attempted 'debug', never does anything.

As far as I can tell, your main problem is mostly caused by the Maple GUI's "document mode". Anyway, I couldn't get it to work. I converted it to "worksheet mode" and 1-D input. I also removed a bunch of unnecessary junk. Both of these were fairly easy to do, and now it works fine.
 

Download Naveshin.mw

The correct expression is
 

`R[1]'`:= ...

Note that the paired outer quotes are different than the solitary inner quote. On my keyboard, the outer quotes are on the upper left, under Esc, and the inner quote is to the left of Enter.

Regarding your followup problem: From a set of sets S you want to extract all L-subsets such that 

1. all pairwise intersections are size 1 (a.k.a., singleton),
2. every element occurs in exactly 2 of the sets.

A procedure for that is

SingletonCliques:= (S::set(set), L::And(posint, Not({1,2})))->
local L2:= L*(L-1)/2, L3:= 2*L-3;
    select(
        s-> local p;
            nops(op~(s)) = L2 and andseq(nops(op~(p)) = L3, p= combinat:-choose(s,2)),
        combinat:-choose(select(nops = L-1, S), L)
    )
:

Again, this is a purely set operation rather than a graph operation, and the edges are atomic with respect to it. You risk causing confusion by mentioning edges or other graph-theoretic aspects.

You wrote "I know this only affects the display," but that's not true. When you see a constant multiplier or divisor distributed over added terms, that's how the expression is actually stored, not just how it displays.

You could use an inert operator to avoid that:

(x^2 - 2*x - 1) %/ 4;

My interpretation of what you want is this:

d:= 2e-4:
theta:= (z,t)->
    -(65.7014900075861*(cos(-4.536529763+45365.29764*z)+.1749541674))
        * exp(-1.603200636*t)
:
T:= [seq](Pi/2+0.1*t, t= 0..7):
plot(
    theta~(z,T), z= 0..d,
    color= COLOR~(HUE, 0.87*T/(max-min)(T)), #0.87: see #1 below
    legend= (t=~ evalf[3](T))
);

#1: 0.87 is a fudge factor to prevent the HUE scale from wrapping 
# around from red at 0 back to red at 1. 

The operation that you asked for is a purely set operation---not a graph operation---and thus the edges are "atomic" with respect to it. In other words, it's irrelevant that the edges are themselves represented as sets; they might as well just be labels. 

This performs your requested operation:

IsPairwiseDisjoint:= (S::set(set))-> local s; evalb(add(nops(s), s= S) = nops(op~(S)))
:
FixedSizeBlocks:= (S::set(set), L::nonnegint)->
    select(IsPairwiseDisjoint, combinat:-choose(S,L))
:

Although that's the operation that you requested, I suspect that there's another set operation that'll actually work better for your graph-theoretic purpose. I'll write that up in a separate entry.

I can't say for sure without seeing it in context, but I suspect that N1 is an initial condition, for example, the value of the formula for n=1. Every recursive formula will have at least one initial condition, although the value of n where it occurs isn't necessarily 1.

The Maple command rsolve can solve some recursion formulas that can be expressed via algebraic equations (in which case they're usually called recurrence relations or difference equatios). Observe here how the specification of initial conditions affects the results:

#No initial condition specified, so Maple chooses N(0) as the initial condition:
rsolve({N(n) = a+N(n-1)}, N(n));
                      N(0) - a + a (n + 1)

#Initial condition specified at n=0:
rsolve({N(n) = a+N(n-1), N(0)=N0}, N(n));
                       N0 - a + a (n + 1)

#Initial condition specified at n=1:
rsolve({N(n) = a+N(n-1), N(1)=N1}, N(n));
                      N1 + a (n + 1) - 2 a

I suspect that the memory explosion occurs when you do the matrix inverse (KL + a0.ML + a1.CL)^(-1). A matrix multiplication of the form X:= A^(-1).B is mathematically equivalent to X:= LinearAlgebra:-LinearSolve(A, B), but the latter is usually much more efficient.

I think that this does what you asked; let me know:

IsLabeledIsomorphicSubgraph:= proc(G1::Graph, G2::Graph)
uses GT= GraphTheory;
local V:= GT:-Vertices(G1), R;
    if {V[]} subset {GT:-Vertices(G2)[]} then 
        R:= GT:-IsSubgraphIsomorphic(G1, GT:-InducedSubgraph(G2, V), 'isomorphism');
        if [R][1] then subs(R[2], GT:-Edges(G1)) else false fi
    fi;
    false
end proc
:

@MaPal93 You asked:

• Are you implying that when Maple's solve() seems to get stuck is because the system has no solution?

No, I'm implying that when there are more equations than unknowns, even just one more, then the system is almost certainly inconsistent, unless you have some strong theoretical reason for believing otherwise. Imagine drawing 3 curves in a plane. It's almost certain that there won't be a single point where all three intersect.

If we can find any inconsistent subset, under any values of parameters, then the whole system is inconsistent. So, the plan below is

1. Let Eq be the set of 12 equations (assumed here to be expressions implicitly equated to 0); let V be the set of 6 variables that you want to solve for.
2. Extract the numerators of the equations; the denominators are irrelevant. Call this EqN.
3. Set any variables in EqN that are not in V to 1, or any other simple, exact value, possibly 0. Call this EqP.
4. Find subsets of EqP of size 6 that contain all variables in V. Call this Eq6.

The above steps are very easy. They'll take 5 seconds or less in total, like this:

EqN:= numer~(Eq):
EqP:= eval(Eq, indets(EqN, name)=~ 1):
Eq6:= select(E-> V subset indets(E, name), combinat:-choose(EqP, 6)):

5. Take any member of Eq6 and try solve on it:

SolE:= {solve}(Eq6[1], V);

6. Remove any solution that contains free variables:

SolE2:= remove(s-> ormap(evalb, s), SolE);

7. If SolE2 is empty, go back to step 5 and pick another set of 6.
8. Take any solution in SolE2 and use it to evaluate the other 6 equations:

EqV:= simplify(eval(EqP minus Eq6[1], SolE2[1]));

9. If EqV = {0}, then, miraculously, the system is consistent, and you have a solution, at least for the parameter valuation used in step 3.
10. If EqV contains anything that is obviously not equal to 0, then the entire system is inconsistent; you're done; period.

• Did you check my script?

Yes. 

I call something like f:= x*cos(x); ff:= x-> f, an attempted indirect parameter reference, the x being the parameter. Those references don't work for procedures created with -> or proc. The most-common way to handle them is ff:= unapply(f, [x]). To be more robust, you should pass x into q.

The quotes that you put around x don't do anything.