The behavior you see is due to using an elementwise operator with multiple container arguments. It's not a bug, and has nothing to do with applyrule. To see what's happening, replace applyrule with foo.

The sublists of the list or Array of "neighbors" should not consist of vertex labels themselves but rather consist of the positive-integer indices into the vertex list that give those labels. So, this is correct:

GraphTheory:-Graph([0,1], [[2], [1]]);
      Graph 2: an undirected unweighted graph with 2 vertices and 1 edge(s)

I'll assume that that is sufficient explanation for you to correct the rest.

The only documentation of this is an extremely poorly written sentence in paragraph 6 of Description at ?GraphTheory,Graph:

• Note that the mapping is an integer mapping that indicates the vertices (if the vertices are labeled as  [1, 2, 3, ...]) or the location of each vertex in the vertex list V.

No example is given.

Also, the word "neighbors" is used incorrectly on this help page when referring to this argument to Graph; it should be "departures". Neighbors are departures plus arrivals. This distinction only matters for directed graphs.

Many graph theory problems can be solved more efficiently by working directly from the adjacency matrix. The following procedure is 20 times faster than your NumberOfK4.

restart:
CountCliques:= proc(G::GRAPHLN, q::posint)
uses GT= GraphTheory, It= Iterator, AT= ArrayTools;
local
    n:= GT:-NumberOfVertices(G),
    A:= AT:-Alias(
        1 -~ rtable(GT:-AdjacencyMatrix(G), datatype= integer[1]),
        [(0..n-1)$2]
    ),
    c, j, ct:= 0
;
    if q>2 then
        for j from 0 to n-1 do A[j,j]:= 0 od;
        for c in It:-Combination(n, q) do
            j:= [seq](c);
            if AT:-IsZero(A[j,j]) then ct++ fi
        od;
        ct
    elif q=1 then n
    else GT:-NumberOfEdges(G)
    fi
end proc
: 
GT:= GraphTheory:
G:= GT:-LineGraph(GT:-CompleteGraph(7)):
CountCliques(G, 3);
                              175

CodeTools:-Usage(CountCliques(G, 4));
memory used=5.73MiB, alloc change=0 bytes, 
cpu time=47.00ms, real time=49.00ms, gc time=0ns

                              105

You have several basic punctuation errors:

1. The last statement in the loop in step 8 ends with "HK." You need to remove the period.
2. You need to put a semicolon after "end do" in the following line.
3. At the end step 13, you need a semicolon after "end do".
4. In step 16, "NFLAG:= 0" should be "NFLAG = 0".
5. At the end of step 18, you need a semicolon after "end if".
6. In step 7, there is an assignment "W3:= W2 + HK*f (T,W2)". You have a blank space after f (it's very hard to see it in the 2D input). You need to remove the space.
7. Same problem with "f (T,W2)" in step 8.
8. Same problem with "f (T,W3)" in step 9.

At this point, after I make the above changes, I am unable to edit or execute the worksheet further. All GUI functions (within this particular worksheet) other than scrolling or saving are blocked; my other open worksheets are not affected. Saving and reopening the worksheet doesn't help.

I don't understand this at all, but I suspect that it has something to do with 2D Input. I've never seen this in the many tens of thousands of hours that I've used Maple over the past 20 years or so. I attach the worksheet here in case someone else can figure out what's going on. I don't think that this has anything to do with your coding; it's a GUI bug.

Download sol_odes.mw

It doesn't work because you don't have write permission to whatever directory (or folder) Maple considers to be your "current" directory. This is likely to be the main directory where the Maple system itself is stored, such as (on Windows) 

currentdir();
                 "C:\Program Files\Maple 2021"

Change that to something that you do have write permission for, such as your desktop. I (personally) would do

currentdir("C:/Users/carlj/desktop"): #Note **forward** slashes (/ not \)!

But you'd need to modify that because you're not "carlj". Once that's done, the save command in your worksheet should work.

It can be done (and, indeed, done surprisingly easily) with the nonlinear programming command Optimization:-NLPSolve like this

restart:
ABS:= (a,b,c)-> abs(a-b*c) <= k*c:
(L,M,U):= seq([cat](x, __, 1..3), x= [l,m,u]);
Cons:= 
    ABS~(
        map(op@index, [L,M,U], [1,1,2]),
        [.67, 2.5, 1.5, 1, 3, 2, 1.5, 3.5, 2.5],
        map(op@index, [U,M,L], [2,3,3])
    )[],
    add(L) + 4*add(M) + add(U) = 6,
    (L <=~ M)[], (M <=~ U)[]
;
Sol:= Optimization:-NLPSolve(k, {Cons}, assume= nonnegative)[2]:
evalf[5](Sol);
 [k = 0.23607,
  l__1 = 0.40747, l__2 = 0.31246, l__3 = 0.16437,
  m__1 = 0.45430, m__2 = 0.36753, m__3 = 0.16437,
  u__1 = 0.54121, u__2 = 0.44972, u__3 = 0.17998]

The Petersen graph has a Hamiltonian path starting at any vertex. Thus a depth-first search for a spanning tree can only find a Hamiltonian path. Therefore, I don't understand why you show a multi-branched tree as your desired outcome of a depth-first spanning tree of the Petersen graph.

[Edit: I've retracted the 2nd sentence in the above paragraph. See my explanation why in a Reply below.]

Anyway, here's a procedure to find spanning trees by depth-first search.

DFSSpanningTree:= proc(G::GRAPHLN, v)
description 
    "Uses depth-first search to construct a spanning tree of "
    "the connected component of G containing v. Edge orientation "
    "and weight (if present) are ignored."
;
uses GT= GraphTheory;
local
    seen:= table(sparse, [v=()]),
    Neighbors:= table(GT:-Vertices(G) =~ GT:-Neighbors(G)),
    stack:= DEQueue(v),
    edges:= rtable(1..0), vertices:= rtable([v]),
    x, y
;
    while not empty(stack) do 
        x:= front(stack);
        for y in Neighbors[x] do
            if seen[y]=0 then
                seen[y]:= (); 
                push_front(stack, y);
                edges,= {x,y}; vertices,= y;
                next 2 #back to `while`
            fi
        od;
        pop_front(stack)
    od;
    GT:-Graph([seq](vertices), {seq}(edges))
end proc
:   
ST:= DFSSpanningTree(
    GraphTheory:-SpecialGraphs:-PetersenGraph(), 1
);
GraphTheory:-DrawGraph(ST);

GraphTheory:-Vertices(ST);
                [1, 2, 3, 4, 5, 8, 7, 6, 10, 9]

The command fnormal can be used to set such small numbers to 0.

Any nonzero scalar multiple of an eigenvector is also an eigenvector for the same eigenvalue.

If we set z = x+I*y for real x, y, then the problem becomes a two-variable optimization problem. Lagrange multipliers is a technique for solving multivariate optimization problems with equality constraints without needing to use the equation to eliminate a variable in the objective (which is often algebraically difficult).

z:= x+I*y:
Obj:= argument(z): #objective function
#Convert constraint equation to equivalent algebraic expression:
Con:= abs(z-6*I) - 3: 
V:= {x,y,m}: #m is the Lagrange multiplier.
crit:= [solve](diff~(evalc(Obj+m*Con), V), V, explicit);
                [ /    1  (1/2)      3  (1/2)      9\   
        crit := [{ m = - 3     , x = - 3     , y = - }, 
                [ \    9             2             2/   

           /      1  (1/2)        3  (1/2)      9\ ]
          { m = - - 3     , x = - - 3     , y = - }]
           \      9               2             2/ ]


eval(z, crit[min[index](eval~(Obj, crit))]);
                         3  (1/2)   9  
                         - 3      + - I
                         2          2  

The code above works just as well with the polar form of complex-to-pair-of-reals correspondence, z = r*exp(I*theta). 

You're looking for a specific numeric solution, so the appropriate command is fsolve, not solve. Yes, there are some similar cases where solve also works, but in those cases you just "got lucky"; fsolve still would've been the better command.

If you want fsolve to restrict itself to E > C without needing to guess an initial value, then replace E=30 with E= C..infinity, and there's no need to use the assume command if you use keyword option continuous as the third argument to int.

Your proposed model function A*f(x,y) (where f(x,y) is the Lowess result from the first pass and A is a real parameter to be fitted) is linear with respect to (wrt) its parameter A, so you can use LinearFit (though I'd recommend adding a constant parameter: A*f(x,y)+B). The fact that this, or any, model is highly nonlinear wrt its independent variables (x and y) is irrelevant. LinearFit should be much quicker than NonlinearFit, and it'll give you the unique values of the parameters that minimize the residuals (in the least-squares sense).

Using the Iterator package, it's probably most efficient to generate directly to a matrix, as in

n:= 4:
<seq(p[], p= Iterator:-CartesianProduct([-1,1]$n))>;

The listlist form can be done by

[seq]([seq](p), p= Iterator:-CartesianProduct([-1,1]$n));

To do the rows in Gray-code order (as in mmcdara's Answer):

2*<seq(p[], p= Iterator:-BinaryGrayCode(n))> -~ 1; #matrix
[seq]([seq](2*p-~1), p= Iterator:-BinaryGrayCode(n)); #listlist

The listlist k is a specification of 12 values of z; the (x,y)-values are determined implicitly, using a scheme similar to what you proposed. The following code will expand any GRID to the (x,y,z)-values of all its points:

ExpandGrid:= proc(
    G:= GRID(range(numeric), range(numeric), listlist(realcons))
)
local a, b, c, d, k:= op(3, G), n:= nops(k)-1, m:= nops(k[1])-1, i, j; 
    (a,b):= (lhs,rhs)(op(1,G)); (c,d):= (lhs,rhs)(op(2,G));
    [seq](
        [seq]([a+j*(b-a)/n, c+i*(d-c)/m, k[j+1][i+1]], j= 0..n), 
        i= 0..m
    )
end proc
:
k:=[[1,2,10],[1,1,8],[0,1,5],[0,0,6]]:
ExpandGrid(GRID(1..2,1..3,k));
        [[           [4      ]  [5      ]           ]  
        [[[1, 1, 1], [-, 1, 1], [-, 1, 0], [2, 1, 0]], 
        [[           [3      ]  [3      ]           ]  

          [           [4      ]  [5      ]           ]  
          [[1, 2, 2], [-, 2, 1], [-, 2, 1], [2, 2, 0]], 
          [           [3      ]  [3      ]           ]  

          [            [4      ]  [5      ]           ]]
          [[1, 3, 10], [-, 3, 8], [-, 3, 5], [2, 3, 6]]]
          [            [3      ]  [3      ]           ]]

This is an Answer to the other question that you asked, the one in the final sentence of your Question, regarding working through a step-by-step solution of a separable ODE from an old worksheet by Harald Pleym.

Maple's int will not automatically perform a definite integral when the integrand contains an unknown function even in cases where the antiderivative seems obvious. Compare:

int(diff(P(t),t), t); #Indefinite integral works.
int(diff(P(t), t), t= 0..T); #Definite integral returns unevaluated.

The reason that it won't do it is that the answer P(T) - P(0) is not necessarily correct^[*1] if P has a discontinuity between 0 and T. Since P is unknown, whether there's a discontinuity is also unknown. You get around this by using option continuous, which tells int to assume that it's continuous (at least between the limits of integration):

int(diff(P(t), t), t= 0..T, continuous);
                          -P(0) + P(T)

The need to use this option may have been introduced after Harald Pleym wrote that worksheet.

If you pass this or any other option to the inert Int, the option will be passed on to int when you use value. You only need to make one change to make the whole worksheet work: In these lines of code:

g=convert(g,parfrac,P(t));
map(x->x*diff(P(t),t),(%));
map(Int,%,t=0..T);

simply change the last line to 

map(Int, %, t= 0..T, continuous);

Footnote: [*1] It's not necessarily correct because different constants of integration may need to be used on the piece containing the lower limit of integration and the piece containing the upper.