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.

Those are called piecewise functions both in standard mathematical English and in Maple. Your phrase "multithreaded function" means something completely different^[*1], and you risk causing great confusion by using that phrase. I recommend that you edit the title of this Question to change "multithreaded" to "piecewise".

Your example function can be specified:

f:= x-> piecewise(x < 0, 1/x, x >= 0, sqrt(x));

See help page ?piecewise.

Footnote: [*1] A multithreaded function is one intended to execute simutaneously in multiple processors sharing the same memory on a single computer. Thus they either do not make assignments to nonlocal memory, or they make those assignments in a carefully controlled way.

The sin(Pi/5) needs to be converted to radical (i.e., algebraic number) form or RootOf form. This does the whole thing, and gives an answer much simpler than Mathematica's:

evala(abs(convert(v, radical)));
                            1  (4/5)
                            - 2     
                            2       

This works also, and may be needed when a "radical" form does not exist:

allvalues(evala(abs(convert(v, RootOf))));

evala is Maple's longstanding "package" of highly sophisticated algebraic-number simplifiers (and other commands). It's not a "package" in the sense of the with command; it's just a collection of procedures. See ?evala.

A `print/...` procedure can be used for any function. This mechanism has been used in Maple for more than 20 years, so it'll certainly work in your version. This doesn't require any typesetting setting, anything from the Typesetting package, any rule assistants, or anything menu-based.

`print/BesselJ`:= (nu,z)-> local J; J[nu](z):

In older versions of Maple, you'll need to use the proc syntax (this is a purely syntactic variation; the functionality is identical):

`print/BesselJ`:= proc(nu,z) local J; J[nu](z) end proc:

This `print/BesselJ` procedure will be called automatically whenever an attempt is made to do 2D output of BesselJ(...).

This would need to be done for every function whose display you want to alter, but it could be done through an initialization file.

Edit: I think that you may be able to fix your issue entirely by using EnableTypesetRule. See, below, my Reply to the Answer by @dharr.