for a to 2 do k||a:= a od;

Here's my encoding of Calvin's algorithm into Maple:

Calvin:= proc(X::procedure, n::posint, gamma::procedure:= (n-> (n+1)^(1/3)))
description "Minimize continuous function with random noise on 0..1";
option
    `Reference: `,
        `https://www.sciencedirect.com/science/article/pii/S0885064X01905746 `,
        `James M. Calvin, "A One-Dimensional Optimization Algorithm and Its `,
        `Convergence Rate under the Wiener Measure", _Journal of Complexity_ `,
        `17, 306-344 (2001)`,
    `Maple author: Carl Love <carl.j.love@gmail.com> 2020-Aug-27`
;     
local 
    T:= Vector(1..2, [0, 1], datatype= float[8]),
    XT:= X~(T), tmin, M, z:= evalf(gamma(0)), k, i, Y, t, tau:= 1
;
    (tmin,M):= `if`(XT[1] < XT[2], [T[1],XT[1]], [T[2],XT[2]])[];
    for k to n do
        Y:= XT -~ (M - z);
        i:= 1 + max[index]( (T[2..] - T[..-2]) /~ Y[2..] /~ Y[..-2]);
        t:= T[i-1] + (T[i]-T[i-1])*Y[i-1]/(Y[i-1]+Y[i]);
        tau:= min(tau, t - T[i-1], T[i] - t);
        ArrayTools:-Insert~([T, XT], i, [t, evalf(X(t))]); 
        if XT[i] < M then (tmin,M):= (t, XT[i]) fi;
        z:= min(z, evalf(gamma(k))*sqrt(tau))
    od;
    [tmin, M]
end proc
:
N:= Statistics:-Sample(Statistics:-RandomVariable(Normal(0, 0.1))):
S:= Array(1..1, datatype= float[8]):
W:= t-> 1 + (t - 0.5)^2 + N(S)[1]:
CodeTools:-Usage(Calvin(W, 350)); 
memory used=25.16MiB, alloc change=0 bytes, cpu time=157.00ms, 
real time=158.00ms, gc time=0ns

             [0.500738492223830, 0.780222456907857]
            

Use

labels= [eta, ` Nu`[a](eta)]

Note that there's a space after the first backquote.

I have only minor changes to suggest, the most significant being replacing the map of verify with a select:

discriminant:= proc(poly::polynom, kv::set(name))
    local f:= poly;    
    while (local vars:= indets(f) minus kv) <> {} do
        f:= (mul@select)(
            t-> kv subset indets(t), 
            op~(1, factors(discrim(f, vars[1]))[2])
        )      
    od;
    f
end proc
:

You use Psi[2] as Psi[2](r*sin(phi)), which makes Psi[2] a function of only one (unnamed) argument. The D(Psi[2]) then refers to the derivative of Psi[2] with respect to that unnamed argument, so it's neither r or phi.

Theoretically, s is a complex parameter, and the inverse Laplace transform is defined via an integration with respect to s over an entire vertical line in the complex plane (so Re(s) is fixed and Im(s) goes from -infinity to infinity). (I'm not saying that this has anything to do with how the inverse transform is computed in practice or in Maple.) So, from this theoretical viewpoint, it makes no sense to assume s<0 or s::real.

See the Wikipedia article "Inverse Laplace transform".

Furthermore, it seems as if you think that assuming s<0 is equivalent to substituting -s for s. It is not. Compare:

sqrt(s^2) assuming s<0;
sqrt((-s)^2) assuming s>0;

Here's a massive simplification of your ApproxSol:

Cw:= <
    .2960214364, .1939557186, 0.2635024425e-1;
    -.2380729574, -0.9481780593e-1, 0.442773709e-1;
    0.7931761947e-1, 0.967928372e-2, -0.2747829289e-1
>:
Ct:= <
    0.7620592980e-10, 0.1532687805e-9, 0.8496500063e-10;
    -0.2362176714e-9, 0.6793447600e-10, 0.2039654777e-9;
    -0.1118492678e-9, -0.5202991920e-10, 0.1128321001e-9
>:
U:= x-> piecewise(x < 0, 0, x < 1, 1):
V:= x-> <1, 4*x-2, 16*x^2 - 16*x + 3>:
ApproxSol:= unapply(w^2 + t + U(w)*U(t)*(Cw.V(w)).(Ct.V(t)), (w,t)):

The program is trivial, but I don't think that you'd want to "display" its results because there are 2^(4^2) = 65,536 of them.

AllRel:= (S::set)-> combinat:-powerset({seq}(seq([i,j], i= S), j= S)):

Any subset of the Cartesian product of a set with itself constitutes a relation on the set.

Using hardware floats and trunc instead of round is much faster:

restart:

f2 := proc(Bins, Probs, N)
  local X, S;
  uses Statistics:
  UseHardwareFloats := false:
  X := RandomVariable(EmpiricalDistribution(Bins, 'probabilities'=Probs)):
  S := Sample(X, N);
  round~(S)
end proc:

fCarl:= (Bins, P::list(realcons), N::posint)->
    (trunc~@[seq]@Statistics:-Sample)(
        Statistics:-RandomVariable(
            'EmpiricalDistribution'(Bins, 'probabilities'= P)
        ),
        N
    )    
:
Probs := [0.10, 0.30, 0.20, 0.25, 0.15]:
Bins  := <50, 60, 70, 80, 90>:

S2 := CodeTools:-Usage(f2(Bins, Probs, 10^5)):
S3:= CodeTools:-Usage(fCarl(Bins, Probs, 10^5)):

Statistics:-Tally~([S||(2..3)]);
memory used=213.87MiB, alloc change=100.00MiB, cpu time=875.00ms, real time=862.00ms, gc time=78.12ms
memory used=8.42MiB, alloc change=1.53MiB, cpu time=47.00ms, real time=34.00ms, gc time=0ns

[[50 = 10011, 60 = 29966, 70 = 20098, 90 = 15035, 80 = 24890], 

  [50 = 10044, 60 = 29912, 70 = 19920, 90 = 15132, 80 = 24992]]

You can't use square brackets in place of parentheses in algebraic expressions. So you need to replace the square brackets that you have in DE1, DE2, and DE3.

I made these three changes:

1. Corrected the spelling of infinity
2. Changed int to Int
3. Changed sum(..., k= 1..1) to eval(..., k= 1).

Now, the lower limit of integration is nonreal and the upper limit is infinity. I don't know what that means mathematically, if the concept is defined at all. And Maple's numerical integration doesn't understand it either.

Overhead: There is a small and fairly constant time cost involved in initializing each thread. This cost is usually called "overhead" (a business analogy). So, to justify starting new tasks, they must have sufficient work to do. The amount of data that you have is much too small for that. When using the Threads:-Task package, it's your responsibility to decide how much work is enough to start a new task. On the other hand, if you use the commands Map, Seq, Add, or Mul from Threads, the command will make that decision for you (and you can help it decide by using the tasksize option).

Globals: Correctly using global variables (your M) with multi-threaded code is very tricky. You should hold off doing that until you're ready for the "advanced class" in shared-memory parallel programming. In many cases, the better solution is to rework the code so that they're not used at all.

Thread safety: The Maple library code had been gradually built up for decades before multi-threaded programming was introduced. The vast majority of that code uses global variables in ways that won't work correctly when multi-threaded. The code that does that is said to be not "threadsafe". See the help page ?index,threadsafe for a list of commands that are known to be threadsafe. You've used the commands floor and sum, which are not threadsafe. You can replace sum with add. You can replace floor with

Floor:= (x::numeric)-> trunc(`if`(x >= 0, x, x-1))

The consequences of using code that's not threadsafe are extremely unpredictable. You can get wrong results without warning, infinite loops, deadlocks ("hangs"), or kernel failures.

You can accommodate parallel edges by adding a degree-2 vertex to them (or to all but one of each group of them). This will not change the count of the number of trails.

For example, in your case, add a vertex 5 to the middle of edge e6, which'll split it into two edges. Likewise add a vertex 6 to the middle of e7.

Add this option to the plotting command: 

axis[2]= [tickmarks= [seq](-1..1, 0.1)]

You might need to reduce the font size to prevent the tick marks from crowding. If so, add another option: 

axesfont= [TIMES, 8]

(for example).

Like this:

seq(seq(U(0,h,m)=0, h= 1..10), m= 0..10)