Apply the expand command to the factored polynomial, the denominator in this case.

I think that you need to have some of Acer's points stated more emphatically: 

1. Your command print(pout) is irrelevant. All that it does is print the filename. 

2. Using the character terminating commands and loops (semicolon or colon) to control the plot output is error prone. You shouldn't use this method.

3. The command that actually produces the final plot that you want saved, display in this case, should be wrapped with a print command. 

1. There's an error in your input of the ODEs: The Pr(...) should be Pr*(...). Maple doesn't catch this type of error when Pr is assigned a numeric value (as it has been); unfortunately, this error simply causes the part in parentheses to be ignored. You have the same problem with alpha(...).

This code is obviously a variation of code originally posted by me. The numerous superfluous parentheses around derivatives indicate that it had been converted to 2D Input and then back to plaintext. I wonder how many hands this code has passed through containing these errors? 

2. You can vary the linestyle using syntax analogous to that used to vary the color. The valid linestyle names are dash, dashdot, dot, longdash, solid, spacedash, spacedot, and ticks. So, for example, linestyle= [solid, dash, dot].

3. This boundary-layer problem uses a finite value as a substitute for an upper boundary at infinity. This should be indicated by a caption on the plots. I think that your 1 is a very bad choice for this value because it obfuscates that the value is a fake upper boundary.

You should not specify a range for eta in the plotting commands. Since this is a BVP, the range will be automatically inherited from the dsolve command.

You should vary the value substituted for infinity to make sure that this does not cause too large a difference in the solutions. The validity of the results is suspect until this has been done. I've read many journal papers on boundary-layer problems that either ignore this point or fail to mention it. Note that the best choice for this value depends on the parameters. This is almost always ignored in these papers.

Maple's mesh-based BVP solvers often become unstable when too large a value is used for a fake infinity in a boundary-layer problem. The solvers can usually detect this instability, but they can't compensate for it; instead they usually return an error message. This is the most difficult issue with solving these boundary-layer problems with Maple.

I'm guessing that the thing that you find most confusing about the procedure is that it doesn't declare its parameters. Instead, it refers to them as args[1], args[2], and args[3]. The args and nargs are keywords: The former is the argument sequence, and the latter is the count of arguments.

Here is my version of the procedure with modern parameter specifications:

BasisGrid:= proc(
    u::{Vector(2, numeric), [numeric, numeric]},
    v::{Vector(2, numeric), [numeric, numeric]},
    ab::{Vector(2, numeric), [numeric, numeric]},
    {pointoptions::list(name, name= anything):=
        [symbol= circle, symbolsize= 20, color= blue]
    }
)
    local combo:= [seq](ab[1]*u+ab[2]*v);
    plots:-display(
        Lamp:-Gridvects([seq](u), [seq](v)),
        Lamp:-Gridvects(
            [1,0], [0,1], max(abs~(combo), 5),
            'vectors'= 'off'
        ),
        plot(combo, 'style'= 'point', pointoptions[]),
        _rest
    )
end proc
:

The _rest is a keyword for arguments that were passed but not declared.

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.