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.

There is no command Remove that would do anything useful with the arguments -10..10, -1. The closest actual command to that is 

remove(`=`, [$-10..10], -1)

which would give you the integers from -10 to 10 with -1 removed. Just using integers won't give you a good plot. But you can do a parametric plot like this:

plot(
    `[]`~(f, g, t=~ [-10..-1.3, -0.7..10]), 
    style= point, symbol= cross, symbolsize= 12, colour= coral
);

Now I'm thinking of the uniform-distribution-of-points problem from your followup.

Use plots:-textplot inside plots:-display to add annotations to a plot. Like this:

InterPlot:= proc(f1,f2)
local 
    F:= [f1,f2], x:= indets(F, And(name, Not(constant)))[], 
    y, v:= [x,y], Int:= solve(y=~ F, v)
;
    plots:-display(
        plot(F, (v=~ (min-1..max+1)~(map(eval~, v, Int)))[]),
        plots:-textplot(
            map2(eval, [v[], v], Int), align= [above,left], color= magenta
        ),
        title= typeset(F)
    )
end proc
:
(print@InterPlot@op)~(
    [x^~[1,3], x^~[2,4], x^~[-1,-2], x^~[2,1/2], [x^3, abs(x^(1/3))]]
):

The next 4 plots are similar and are omitted.