When the while condition is evaluated the first time, n is not assigned a value.  So the condition n < nmax cannot be evaluated to true or false. Note that the call to evalb is not needed.

Duncan's suggestion to use the CodeTools package is a good one.  To use these tools, the code needs to be in form of executable procedures, rather than top-level commands.  If the current code is one large worksheet/document, you might find it convenient to export it as Maple input, then wrap the contents of that text file into a Maple procedure.  Depending on how the orginal worksheet was coded, the top-level variables may have to be declared as globals, but if you can avoid doing so it may run faster since global variabled are fully evaluated each time they are accessed. 

As an example, suppose you started with a simple worksheet that looked like

L := NULL:  # this is a lousy way to create the sequence 1,4,9,16,...,100
for i to 10 do
    L := L, i^2;
end do;

After exporting to a text file and wrapping in a procedure you might have a file named createL.mpl with the content

CreateL := proc()
local i, L;       # declare these variables as locals
   L := NULL;
   for i to 10 do
      L := L, i^2;
   end do;
   return L; 
end proc:

Inside of Maple you would then do

read "createL.mpl":
with(CodeTools:-Profiling):
Profile(CreateL):
L := CreateL();
                   L := 1, 4, 9, 16, 25, 36, 49, 64, 81, 100

PrintProfiles();
CreateL
CreateL := proc()
local i, L;
     |Calls Seconds  Words|
PROC |    1   0.000    174|
   1 |    1   0.000      0| L := NULL;
   2 |    1   0.000      0| for i to 10 do
   3 |   10   0.000    174|   L := L, i^2
                            end do;
   4 |    1   0.000      0| return L
end proc

stackmatrix and augment are part of the deprecated linalg package, they do not work with rtable based matrices (created with Matrix). They aren't needed with rtables.  You can do, for example,

(**) A := Matrix(2, symbol='a');
                                                             [a[1, 1]    a[1, 2]]
                                                        A := [                  ]
                                                             [a[2, 1]    a[2, 2]]

(**) B := Matrix(2, symbol='b');
                                                             [b[1, 1]    b[1, 2]]
                                                        B := [                  ]
                                                             [b[2, 1]    b[2, 2]]

(**) <A|B>;                     
                                               [a[1, 1]    a[1, 2]    b[1, 1]    b[1, 2]]
                                               [                                        ]
                                               [a[2, 1]    a[2, 2]    b[2, 1]    b[2, 2]]

(**) <A,B>;                     
                                                          [a[1, 1]    a[1, 2]]
                                                          [                  ]
                                                          [a[2, 1]    a[2, 2]]
                                                          [                  ]
                                                          [b[1, 1]    b[1, 2]]
                                                          [                  ]
                                                          [b[2, 1]    b[2, 2]]

What's the problem?

As an example, suppose you are looking for all variables of the form x(t), y(t), etc.  One way to do that with ?indets is

(**) ex := diff(x(t),t,t) + sin(y(t)) + cos(t);
                                             / 2      \
                                             |d       |
                                       ex := |--- x(t)| + sin(y(t)) + cos(t)
                                             |  2     |
                                             \dt      /

(**) indets(ex, 'And(anyfunc(identical(t)),Not(known))');
                                                   {x(t), y(t)}

Note the use of the structured type (see ?type,structured) to select the appropriate variables.  The type known, alas, is not documented, but you can see its assignment by doing

(**) `type/known`;
                                      And(function, PDEtools/known_function)

The purpose of it here is to eliminate cos(t), which is a known function and so generally not suitable as a variable.

How about using ?LinearAlgebra[JordanForm]? Hmm. I wasn't paying attention, you only want to permute columns...

The limit depends on the sign of T.  Try

limit(exp(x/T),x=infinity) assuming T>0;

I don't believe there is a way to use more than one core; dsolve/numeric is not designed to be multi-threaded.

Have you tried using a different solver, say rosenbrock? 

Pagan's suggestion is probably what you want.  However, if you happen to be sorting a large list, a faster approach is to sort with attributes. For example,

 map(attributes,sort([seq(setattribute(evalf(x),x),x in a)]));
                          1/2   1/2     1/2    1/2    1/2
                       [-5   , 2   , 2 2   , 26   , 39   ]

Here is a simple procedure that converts a Vector to a tilde matrix:

tilde := proc(V :: Vector)
local T;
    T := Matrix(3, 'shape=skewsymmetric'):
    T[3,2] := V[1];
    T[1,3] := V[2];
    T[2,1] := V[3];
    return T;
end proc:

V := <x,y,z>:
tilde(V);
                                            [0     -z    y ]
                                            [              ]
                                            [z     0     -x]
                                            [              ]
                                            [-y    x     0 ]

Because H(z) is not defined, I don't see how you can get a numerical solution.

If that is all you want to do, a fairly simple way is to the following one-liner:

plots:-display(seq(plot(2*x+i/2,x),i=0..10));

To do this with a loop, you should save each plot in a table and then make a call to display to plot them at once.

for i from 0 to 10 do
     y := i/2;
     f := 2*x + y;
    P[i] :=  plot( f, x, color=black);
end do:
plots:-display(seq(P[i], i=0 .. 10));

The  custom component generates one positive pulse of a sinewave with frequency 200 and amplitude 235.  One way to generate the same signal without a custom component is to multiply the output of a sinewave with a step function that transitions from 1 to 0 after a half period of the sine wave.

The most likely problem is that the number of lines in the file is not what you expect. Change your code to

N := iquo(nops(Data),3);
Data := convert(Data,Array);

Also, rather than using the deprecated array structures, and later converting to Vectors, just initially assign them as Vectors:

coords_x := Vector(1..N):
coords_y := Vector(1..N):
Ex := Vector(1..N):
Ey := Vector(1..N):
Ez := Vector(1..N):

and remove the final conversions.  Here's the code I used:

restart:
with(plots):
file := fopen("900MHz.dat",READ):
Data := readdata(file,7):
N := iquo(nops(Data),3);
Data := convert(Data,Array);

coords_x := Vector(1..N):
coords_y := Vector(1..N):
Ex := Vector(1..N):
Ey := Vector(1..N):
Ez := Vector(1..N):
t := 1:
for i to N do
    coords_x[i] := Data[t,1]:
    coords_y[i] := Data[t,2]:
    Ex[i] := Data[t,6]:
    Ey[i] := Data[t + 1,6]:
    Ez[i] := Data[t + 2,6]:
    t := t + 3:
end do:
field := Vector(1..N):
for i from 1 to N do
    field[i] := sqrt(Ex[i]^2 + Ey[i]^2 + Ez[i]^2):
end do:
pointplot3d(coords_x, coords_y, field, 'axes=normal, title="Near field plot"');

Trying simplifying the problem and then comparing with a known result.  For example, remove the losses from the end, so the boundary conditions are 

IBC := {NULL
        , -lambda*(D[1](T))(L, t) = 0
        , T(0, t) = 1573
        , T(x, 0) = T0
       };

The time-constant of the model is rho*cp/lambda*L^2 ~ 12 hours.  A distributed RC model should reach about 90% of the input step after one time-constant.  Here pdsolve gives 88.6%.  So that seems to be working.  You should check that against your Matlab/FEA simulations.