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)

@Bohdan Unless I've made some error in my simplification of the coding of your system (which is entirely possible), the system that I get is full rank, so the Moore-Penrose and regular inverses are the same. This is my recoding, please check it:

restart:
local gamma: #Overrides Maple's predefined gamma constant.
EQ:= P[j,r,u] = Sum(
    pi[i,j,r,u]*(gamma[i,r]*rho[i,r] + Sum(gamma[i,s,r]*P[i,s,r], s= 1..S)),
    i= 1..J
);
eq:= value(eval(EQ, [J= 3, S= 2])):
sys:= [seq](seq(seq(eq, j= 1..3), r= 1..2), u= 1..2);
Ps:= indets(sys, specindex(P));

I strongly suspect that the "hard time" or sluggishness that you're having further processing the output has nothing to do with its mathematical complexity but rather with the poor ability of Maple's Standard GUI (which is essentially a separate program (written in Java) from Maple's mathematical components, and runs as a separate process) to format that output for display on-screen. There are several things that you can do to ameliorate this. You can do any subset of these; they are independent:

1. Do not copy and paste output. Do not use equation labels. Instead, in the input command that generates the output, assign the output to a variable, like I show above. Strive to make every command directly typed rather than generated from prior output.
2. If you don't need to see the output of a command, end it with a colon. 
3. If you've seen some output but no longer need to see it, use Ctrl-Delete to delete it. Since it's all stored in variables (by point #1), it can be redisplayed without recomputation whenever you need to see it again.
4. Do not select large amounts of output with the mouse. Thus, don't use context menus on them.
5. In extreme cases, give the command interface(prettyprint= 0): before the command that generates problematic output. You can reset it to normal by interface(prettyprint= 3):