By avoiding n^2 list building and using a few other refinements, I reduced the run time of Routes from about 20 seconds to 1/2 second. Here is the modified code

Routes:=proc(N::posint, n::nonnegint)

global T;
local i,j,L, Rule;

    if n>=N^2 then
        T:=[]:
        0;
    else
        Rule := proc(K)  # Continuation of the route by 1 step
        local S, k, r, rk, p, pts, j;

            S := table();
            j := 0;
            k := nops(K[1]);

            for r in K do
                # Assign pts the points to consider
                if r[k]=[1, 1] then
                    # Bottom left corner
                    pts := [[1, 2], [2, 2], [2, 1]];
                elif r[k]=[1, N] then
                    # Top left corner
                    pts := [[1, N-1], [2, N-1], [2, N]];
                elif r[k]=[N, N] then
                    # Top right corner
                    pts := [[N-1, N], [N-1, N-1], [N, N-1]];
                elif r[k]=[N, 1] then
                    # Bottom right corner
                    pts := [[N-1, 1], [N-1, 2], [N, 2]];
                elif r[k,1]=1 and r[k,2]<>1 and r[k,2]<>N then
                    # Left side
                    pts := [[1, r[k,2]-1], [2, r[k,2]-1], [2, r[k,2]-1], [2, r[k,2]+1], [1, r[k,2]+1]];
                elif r[k,2]=N and r[k,1]<>1 and r[k,1]<>N then
                    # Top side
                    pts := [[r[k,1]-1, N], [r[k,1]-1, N-1], [r[k,1], N-1], [r[k,1]+1, N-1], [r[k,1]+1, N]];
                elif r[k,1]=N and r[k,2]<>1 and r[k,2]<>N then
                    # Right side
                    pts := [[N, r[k,2]+1], [N-1, r[k,2]+1], [N-1, r[k,2]], [N-1, r[k,2]-1], [N, r[k,2]-1]];
                elif r[k,2]=1 and r[k,1]<>1 and r[k,1]<>N then
                    # Bottom side
                    pts := [[r[k,1]-1, 1], [r[k,1]-1, 2], [r[k,1], 2], [r[k,1]+1, 2], [r[k,1]+1, 1]];
                elif r[k,1]<>1 and r[k,1]<>N and r[k,2]<>1 and r[k,2]<>N then
                    # Inside
                    pts := [[r[k,1]-1, r[k,2]-1], [r[k,1]-1, r[k,2]], [r[k,1]-1, r[k,2]+1], [r[k,1], r[k,2]+1], [r[k,1]+1, r[k,2]+1], [r[k,1]+1, r[k,2]], [r[k,1]+1, r[k,2]-1], [r[k,1], r[k,2]-1]];
                fi;
                rk := r[1..k-1];
                j := j+1;
                S[j] := seq(`if`(member(p,rk)
                                 , NULL
                                 , [op(r),p]
                                ), p = pts);
            od;
            [seq(S[j], j=1..j)];
        end proc;

        L:=[seq(seq([[i, j]], j=1..N), i=1..N)];
        T:=(Rule@@n)(L);  # List of all the routes
        nops(T);  # Number of all the routes
    fi;
end proc:

m := 5:
LinearAlgebra:-BandMatrix([[seq(a(k),k=2..m-1)],[seq(b(k),k=1..m-1)],[seq(a(k),k=1..m-1)]]);
                                      [b(1)    a(1)     0       0       0  ]
                                      [                                    ]
                                      [a(2)    b(2)    a(2)     0       0  ]
                                      [                                    ]
                                      [ 0      a(3)    b(3)    a(3)     0  ]
                                      [                                    ]
                                      [ 0       0      a(4)    b(4)    a(4)]

with(Statistics):
X := RandomVariable(Binomial(10,1/2));
ProbabilityFunction(X,3);
                                      15
                                      ---
                                      128

The problem lies in the first ln: (1+(XL)^(2))(1+(YL)^(2).  There should be a multiplicative operator between the two additive terms. Without that, Maple interprets this as function application.  I'll make a change to AlgEquation so that it returns the offending terms when it detects an error.

Edit -> Remove Output (pick either Worksheet or Selected Region)

You can use the Iterator package to solve this puzzle with nested parentheses.  To step through all possible groupings,  use the BinaryTrees iterator.  The following appliable module does that. Understanding precisely how it does it might take a bit of study. This isn't ideal; one really should prune the trees to prevent, for example, both the following branches (1 + 2) + 3 and 1 + (2 + 3) from appearing.  Doing so isn't hard but I'll leave that as an exercise.  Looping through all possibilities takes a while; here I illustrate the operation using n=6, which limits the tree to six internal nodes (7 leaves), so the digits are from 1 to 7 rather than 1 to 9.

DigitalPuzzle := module()
export ModuleApply;
local Expr, Format;

    # Given the L and R Arrays from the BinaryTree iterator,
    # construct an expression corresponding to the tree,
    # where o[v[i]] is the arithmetic operator of the i-th
    # internal node (the v-Array is used to change the operators
    # for a give tree).
    Expr := proc(L,R)
    local leaf,prefix;
    global o,v;
        leaf := 0;
        prefix := proc(i)
            if i=0 then leaf := leaf+1;
            else 'o'['v'[i]](prefix(L[i]),prefix(R[i]));
            end if;
        end proc;
        prefix(1);
    end proc:

    # Given the L and R arrays of the current tree, and an array, v,
    # that maps the internal nodes to the arithmetic operators (stored
    # in o), return a string corresponding to the desired equality.
    # This is only used for formatting a correct result; speed is
    # not a significant issue.
    Format := proc(L,R,v,o,targ)
    local leaf,infix;
        leaf := 0;
        infix := proc(i)
            if i=0 then leaf := leaf+1;
            else sprintf("(%A %A %A)", infix(L[i]), o[v[i]], infix(R[i]));
            end if;
        end proc;
        sprintf("%s = %a", infix(1), targ);
    end proc:

    ModuleApply := proc(n::posint, targ:=100)
    local A,Accept,cnt,ops,Op,BT,LR,s;
    uses Iterator;

        # Array of arithmetic operators
        ops := Array(0..3,[`+`,`*`,`-`,`/`]);

        # Template for the accept predicate.  The try/catch is
        # needed to handle division by zero.
        Accept := proc(v,o:=ops)
            try
                evalb(_ex=targ);
            catch:
                false;
            end try;
        end proc;

        A := ()->NULL; # dummy procedure that is replaced

        # Construct an iterator that generates all possible values of
        # arithmetic operators (as an Array with values from 0 to 3
        # corresponding to the four operations), but accepts (outputs)
        # only those that meet the criteria.
        Op := MixedRadixTuples([4$n], 'accept'=A);

        # Construct an iterator that generates all binary trees
        # with n-internal nodes (and n+1 leaves).
        BT := BinaryTrees(n);

        cnt := 0; # success counter

        # Loop through all binary trees
        for LR in BT do
            # Reassign A, which is the accept predicate in
            # the Op[erator] iterator, specializing it to the
            # selected tree.
            A := subs(_ex=Expr(LR),op(Accept));

            # Loop through all possibilities of assigning
            # the arithmentic operators to the internal nodes of the
            # tree, keeping only those that evaluate to the target.
            for s in Op do
                cnt := cnt+1;
                printf("%s\n", Format(LR,s,ops,targ));
            end do;
            reset(Op);
        end do;
        cnt;
    end proc:
end module:

CodeTools:-Usage(DigitalPuzzle(6));
(1 + (((((2 + 3) * 4) * 5) + 6) - 7)) = 100
((1 + ((((2 + 3) * 4) * 5) + 6)) - 7) = 100
(((1 + (((2 + 3) * 4) * 5)) + 6) - 7) = 100
(((1 - 2) + 3) + (((4 * 5) - 6) * 7)) = 100
(1 - ((2 * 3) - (((4 + 5) + 6) * 7))) = 100
(1 - ((2 - 3) - (((4 * 5) - 6) * 7))) = 100
((1 - (2 * 3)) + (((4 + 5) + 6) * 7)) = 100
((1 - (2 - 3)) + (((4 * 5) - 6) * 7)) = 100
((1 - 2) + (3 + (((4 * 5) - 6) * 7))) = 100
((1 * 2) * (3 + (((4 + 5) * 6) - 7))) = 100
(1 * (2 * (3 + (((4 + 5) * 6) - 7)))) = 100
(1 - (2 - (3 + (((4 * 5) - 6) * 7)))) = 100
((1 * 2) * ((3 + ((4 + 5) * 6)) - 7)) = 100
(1 * (2 * ((3 + ((4 + 5) * 6)) - 7))) = 100
((1 + ((2 * 3) * 4)) * ((5 + 6) - 7)) = 100
((1 + (2 * (3 * 4))) * ((5 + 6) - 7)) = 100
(1 - ((2 * 3) - ((4 + (5 + 6)) * 7))) = 100
((1 - (2 * 3)) + ((4 + (5 + 6)) * 7)) = 100
(1 + ((((2 + 3) * (4 * 5)) + 6) - 7)) = 100
((1 + (((2 + 3) * (4 * 5)) + 6)) - 7) = 100
(1 + ((((2 + 3) * 4) * 5) + (6 - 7))) = 100
((1 + (((2 + 3) * 4) * 5)) + (6 - 7)) = 100
(((1 + ((2 + 3) * (4 * 5))) + 6) - 7) = 100
((1 - (2 * 3)) * ((4 * 5) * (6 - 7))) = 100
((1 - (2 * 3)) * ((4 * 5) / (6 - 7))) = 100
((((1 - (2 * 3)) * 4) * 5) * (6 - 7)) = 100
((((1 - (2 * 3)) * 4) * 5) / (6 - 7)) = 100
(1 + (((2 + 3) * (4 * 5)) + (6 - 7))) = 100
((1 + ((2 + 3) * (4 * 5))) + (6 - 7)) = 100
((1 + ((2 * 3) * 4)) * (5 + (6 - 7))) = 100
(((1 - (2 * 3)) * (4 * 5)) * (6 - 7)) = 100
(((1 - (2 * 3)) * (4 * 5)) / (6 - 7)) = 100
(((1 - (2 * 3)) * 4) * (5 * (6 - 7))) = 100
(((1 - (2 * 3)) * 4) * (5 / (6 - 7))) = 100
((1 + (2 * (3 * 4))) * (5 + (6 - 7))) = 100
((1 - (2 * 3)) * (4 * (5 * (6 - 7)))) = 100
((1 - (2 * 3)) * (4 * (5 / (6 - 7)))) = 100
memory used=370.72MiB, alloc change=24.00MiB, cpu time=6.43s, real time=6.60s
                                                                    37

@JotaTR You should contact Maple Support. Did you install the MapleSim 5.02 update?  That will probably fix this.

This is known as Dudney's Digital Century puzzle.  Coincidentally, I included a simplified version of it as an example in the Iterator package I added to the Maple Application Center a few weeks ago.  Actually, that isn't quite accurate; the example is in an update to the package that I haven't yet uploaded.  I've been adding various features to it over the holidays.  Here is how I solved the simplified version (restricting the operators to just  addition and  multiplication, with no parentheses).

with(Iterator):

# Assign a procedure that converts a boolean Vector to the relevant
# string.

VecToStr := proc(V)
local i;
    cat("",seq('(i,`if`(V[i]=0,"+","*"))',i=1..8),"9=100")
end proc:

# Construct an iterator that steps through all 8-bit boolean Vectors,
# keeps those for which the equation is true, and returns the
# transformed string.
iter := BinaryGrayCode(8
                       , 'accept' = (V -> evalb(parse(VecToStr(V))))
                       , 'transformer' = VecToStr
                      ):
# Use the iterator to generate all the solutions.
seq(s, s=iter);
   "1*2*3*4+5+6+7*8+9=100", "1*2*3+4+5+6+7+8*9=100", "1+2+3+4+5+6+7+8*9=100"

Here's a quick extension to handle all four arithmetic operations, but without parentheses.

VecToStr := proc(V)
local i,ops;
    cat("",seq('(i,"+*/-"[V[i]+1])',i=1..8),"9=100")
end proc:

iter := MixedRadixTuples([4$8]
                         , 'accept' = (V -> evalb(parse(VecToStr(V))))
                         , 'transformer' = VecToStr
                        ):
seq(s, s=iter);
"1+2+3+4+5+6+7+8*9=100", "1+2+3-4*5+6*7+8*9=100", "1+2*3+4*5-6+7+8*9=100",
    "1+2*3*4*5/6+7+8*9=100", "1+2-3*4+5*6+7+8*9=100", "1+2-3*4-5+6*7+8*9=100",
    "1*2*3+4+5+6+7+8*9=100", "1*2*3*4+5+6+7*8+9=100", "1*2*3*4+5+6-7+8*9=100",
    "1*2*3-4*5+6*7+8*9=100", "1-2+3*4*5+6*7+8-9=100", "1-2+3*4*5-6+7*8-9=100",
    "1-2*3+4*5+6+7+8*9=100", "1-2*3-4+5*6+7+8*9=100", "1-2*3-4-5+6*7+8*9=100"

The structure you have created is a ?list, not an ?Array. Use ?numelems or ?nops to compute its length.

The conversion is correct but is for the relative temperature. Use ?convert/temperature to convert absolute temperatures.

_params can only be indexed by one element (_params['b','c'] makes no sense).  Changing _params['b','c'] to _params['b'] fixes your procedure.

There are others ways to achieve the same thing.  Here are two

f1 := proc() local x; add(x, x = args) end proc: 

f2 := proc(a:=0, b:=0, c:=0) a+b+c end proc:

The key is that you are using ?PDEtools[declare], which uses Maple's ?print facility to modify the display of the output. Fortunately, there is a simple hack to allow you to access this stuff directly.

 sprintf("%A & %A", op(eval(S,[diff=`print/diff`,`declare/display`[2][]])));

That should give you what you want. Note, however, that, in the GUI, the ampersand is displayed by itself (without the "amp;").  I'm guessing that isn't really a problem.  The string actually consists of what you want; to verify that, do lprint(%).  I'm thinking you want to send this to some other application, so that should suffice.

Addendum

I suppose some explanation would be in order.  The ?eval command is used to replace diff(x(t),t) with `x'` and x(t) with x.   The former replacement (diff(t),t) --> `x`')  is handled by replacing the diff function with `print/diff`, which returns the form used when displaying the result.  The latter replacement is handled by using the list of replacements, stored in `declare/display`[2].  It is used by PDEtools[declare] to build the replacements.

The %A format string is used in ?sprintf, rather than %a, to elide the backquotes around `x'`.

Generalization

The technique described above works, but requires some knowledge of PDEtools[declare] to figure out. Actually, my knowlege of PDEtools[declare] is minimal; I was able to quickly figure out the hackery above by using a Maple debugger. In hindsight, that wasn't necessary. That is, it is possible to achieve the same goal in a more general way.  Here is a simple procedure that replaces Maple subexpressions that with their printed equivalent.

toPrint := proc(ex)
local n,p;
    p := select(type, [anames(procedure)], 'suffixed(`print/`)');
    n := map(substring, p, 7..-1);
    eval(ex, n =~ p);
end proc:

It can be used in the above by doing

sprintf("%A & %A", op(toPrint(S)));
         "x' = y & y' = x^2"

To sort the columns, without creating a permutation Matrix, you could do

(m,n) := upperbound(M):
for j to n do
    M[..,j] := sort(M[..,j]);
end do:

Here's another approach for generating permutation lists/vectors.

There is a syntax error in the  large expression. I'm not sure how you want to resolve it.  The exp((.28*...   near the end is never completely closed. Removing one of the open parentheses allows Maple to solve it very quickly.

A usual method for solving this sort of problem is to use Lagrange multipliers. That avoids having to solve the constraint for one of variables. Maple provides a routine for this,

(**) V := x*y*z:                                                          
(**) L:=4*x+2*y+2*z:                                                      
(**) Student[MultivariateCalculus][LagrangeMultipliers](L,[V-32],[x,y,z]);
                  1/2              1/2            1/2
[2, 4, 4], [-1 + 3    I, -2 + 2 I 3   , -2 + 2 I 3   ],

           1/2              1/2            1/2
    [-1 - 3    I, -2 - 2 I 3   , -2 - 2 I 3   ]