Remove the option y= 0..100.

ODE2:=(Diff(T(t), t) = 1/59*(18-T(t))):
DEtools[DEplot]( ODE2, T(t), t=-4..4, {[T(0)=88]});

Your problems are caused by your combining several with commands into one statement. You have

with(plots), with(ColorTools), with(LinearAlgebra), with(RandomTools), with(ExcelTools);

You should change that to

with(plots); with(ColorTools); with(LinearAlgebra); with(RandomTools); with(ExcelTools);

The results of chaining with statements seems unpredictable. Note the following proviso from ?with :

The with command is effective only at the top level, and intended primarily for interactive use. Because with operates by using lexical scoping, it does not work inside the bodies of procedures, module definitions, or within statements. [emphasis added]

Although I can't tell what's wrong with your code until you upload it, here's an example of correctly doing what you want to do. Maybe you'll be able to correct your code from this.

restart:
with(VectorCalculus):
SetCoordinates(cartesian[x,y]):
V:= VectorField(< x+y, x-y >):
F:= (X,Y)-> VectorCalculus:-Norm(VectorCalculus:-evalVF(V, < X,Y >)):
F(1,2);

Here's how to apply the weights from your main graphs to your subgraphs constructed from the paths.

ApplyWeights:= proc(G1::GraphTheory:-Graph, G2::GraphTheory:-Graph)
#Apply the weights from G1 to G2.
uses GT= GraphTheory;
local W:= GT:-WeightMatrix(G1, copy= false);
     GT:-Graph(
          map(e-> [e, W[e[]]], GT:-Edges(G2, weights= false))
     )
end proc:          

G4:= ApplyWeights(G1, G3);

G5:= ApplyWeights(G2, G3);

Here's how to do it. I used essentially the same example as I did for your previous question, except that I made it a directed graph instead of undirected, because I think that's more relevant to your application.

Download Two_Paths.mw

You have a line

F[total] = F[C4H10](W) + ...

which should be changed to

F[total]:= F[C4H10](W) + ....

After that, you cannot use F[total](0) as a variable. So change the three occurrences of F[total](0) to something else, such as F_total_0.

I made these changes, and the dsolve and plots come out perfectly.

Here's an example based on an example found on the help page ?DynamicSystems,TransferFunction .

ss_a:= Matrix([[1, 2], [0, 4]]):
ss_b:= Matrix([[3, 7], [9, 6]]):
ss_c:= Matrix([[5, 6], [5, 2]]):
ss_d:= Matrix([[0, 0], [0, 0]]):
sys:= DynamicSystems:-TransferFunction(
     ss_a, ss_b, ss_c, ss_d,
     discrete, sampletime=0.001, systemname="MIMO system"
);

sys:-tf[1,1];

sys:-tf[1,2];

And likewise for sys:-tf[2,1] and sys:-tf[2,2]. In order for me to give you more details, you'll need to post the code that led to the RTABLEs you showed.

The command to find the minimal-weight path (assuming the weights are nonnegative) is GraphTheory:-DijkstrasAlgorithm, not GraphTheory:-ShortestPath.

Here is a complete solution to your Questions 1 and 2:

1. Construction of the random weighted graph:

restart:
randomize(3): #Needed to get a nontrivial shortest path.
N:= 100: #Number of vertices.
E:= 2300: #Number of randomly chosen edges.
Edges:= map(`{}`@op, combinat:-randcomb(combinat:-choose(N,2), E)):
Weights:= RandomTools:-Generate(list(float(method= uniform), E)):
WeightedEdges:= {zip(`[]`, Edges, Weights)[]}:
G:= GraphTheory:-Graph(WeightedEdges);

2. Find the minimal-cost path and its cost:

Sol:= GraphTheory:-DijkstrasAlgorithm(G, 1, N);

(Path,Cost):= Sol[];

Path:= [seq({Path[k], Path[k+1]}, k= 1..nops(Path)-1)]:

3. Find the maximal-cost edge on that path:

W:= GraphTheory:-WeightMatrix(G, copy= false):
MaxWeight:= -infinity:
for e in Path do
     if W[e[]] > MaxWeight then
          MaxWeight:= W[e[]];
          MaxEdge:= e
     end if
end do:

4. Construct the subgraph that has that edge removed:

G1:= GraphTheory:-Graph(WeightedEdges minus {[MaxEdge, MaxWeight]});

Here's how to solve your followup request---using larger initial sets. This will handle the case in the first Answer also.

#Step 0:
B1:= A1+A2:  B2:= A1+A3:  B3:= A3+A4:  B4:= A2+A4:
C1:= `+`(A||(1..4)):
R:= {A||(1..4), B||(1..4), C1}:

#Step 1:

The procedure X3 takes any set and produces a set with three times as many elements, using the process that you described.

X3:= proc(R::set)
local
     Names:= indets(R,name),
     S0:= [R[]],
     A,
     S01:= subs(seq(A= cat(A,1), A in Names), S0),
     S02:= subs(seq(A= cat(A,2), A in Names), S0)
;
     {S01[], S02[], zip(`+`, S01, S02)[]}
end proc:

S1:= X3(R):  T1:= X3(S1):  #etc.
nops(T1); #Confirmation
     81

#Step 2:

This is exactly the same as in the first answer, with S1 replaced with T1, etc. Is that enough detail for you to finish it?

Here's how to do it. I've coded it all as top-level code rather than as a procedure that produces each iteration. I removed all the underscores from your notation just to spare myself from typing them.

restart:

#Step 0:
B1:= A1+A2:  B2:= A1+A3:  B3:= A3+A4:  B4:= A2+A4:
C1:= `+`(A||(1..4)):
R:= [A||(1..4), B||(1..4), C1]:

#Step 1:
S01:= subs(seq(A||i= A||i||1, i= 1..4), R):
S02:= subs(seq(A||i= A||i||2, i= 1..4), R):
S1:= {S01[], S02[], zip(`+`, S01, S02)[]}:
nops(S1); #Confirmation
     27

#Step 2:
z:= P-> `+`(P[]) mod 2:
S2:= map(z, combinat:-choose(S1,2)) minus S1:
nops(S2); #Confirmation
     162

I adapted a cartesian-product procedure from the MaplePrimes archives:

CartProdSeq:= proc(L::seq(set))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
    eval({subs(Seq=seq, foldl(Seq
                              , [cat(i,1..nargs)]
                              , seq(cat(i,j)=L[j],j=nargs..1,-1)
                             ))});
end proc:

And with that machinery built, your main loop is simply 

Union:= S1:
for n from 3 while S||(n-1) <> {} do
     Union:= Union union S||(n-1);
     S||n:= map(z, CartProdSeq(S1,S||(n-1))) minus Union
end do;

I omitted the lengthy output. But S3 has 66 elements and S4 is empty.

surfdata requires that an Array of [x,y,z] data be 3D.  Since your M is a 441 x 3 Array, and since 441 = 21^2, I assumed that M corresponded to a 21 x 21 MESH.

M:= ExcelTools:-Import("C:/Users/Carl/desktop/Mdata.xlsx"):
M3:= ArrayTools:-Alias(M, [1..21, 1..21, 1..3]):
plots:-surfdata(M3);

If you evaluate r(X) at randomly chosen numeric values of X, then the value (i.e., the rank) will "almost always" be maximal. To quantify "almost alawys", use random integer X and compute the rank over a large finite field, where the probability that the rank of the numeric matrix is not maximal is a very small but easily computable number. Repeat independent trials until the product of the probabilities is below any epsilon that you want. Use enough different fields to lift the results to the reals.

Make f this:

f:= proc(t,M)
local n;
     if not [args]::list(numeric) then 'procname'(args)
     else evalhf(add(-4/Pi/(2*n+1)*sin((2*n+1)*t), n= 0..M))
     end if
end proc:

This involves passing M to f in addition to t. Then the plot call looks this:

plot(f(t,100), t= 0..2*Pi);

This runs in less than 1/10 second, using M=100.

Issue the command

Units:-UseSystem('CGS');

Names are by default symbols in Maple; you don't need to do anything to make them symbols. They stay symbols until you assign something to them. Your problem is caused by a lack of assumptions on your constants.

You do not need to enter with(Statistics) twice. You should not modify the global sequence constants. That sequence is for constants that have a definite value, like Pi.

restart:
with(Statistics):
h:= RandomVariable(Exponential(H)):
g:= RandomVariable(Exponential(G)):
assume(Pr > 0, Ps > 0):
PDF(Ps*g*h/(Pr*g+2), t);