This is a known bug of the Standard GUI that has been reported since Maple 13. See e.g. this thread and the linked one.

It happens that subs is not mapping over M:

> map2(subs,ans,M);
                                [%1    %2    %3]
                                [              ]
                                [%1    %2    %3]
                                [              ]
                                [%1    %2    %3]
                                [              ]
                                [%1    %2    %3]
                                          1
                          %1 := ----------------------
                                1 + exp(-0.1603426501)
                                          1
                           %2 := --------------------
                                 1 + exp(1.734601055)
                                          1
                          %3 := ----------------------
                                1 + exp(-0.7537718024)

See this thread and those linked therein.

What is needed is a good, flexible enough complete square command at top level. One of the problems with Student:-Precalculus:-CompleteSquare is its restiction on the type of the arguments after the first one (basically names or functions). So, a workaround is writing a wrapper procedure that handles more general expressions for those arguments, like in the following sketch:

> completesq:=proc(ex::algebraic,l::{algebraic,list(algebraic)})
> local l1:=freeze~(l),
> s:=zip(`=`,l,l1),
> s1:=foldr(algsubs,expand(ex),`if`(s::list,op(s),s)),
> s2:=Student:-Precalculus:-CompleteSquare(s1,l1);
> thaw(s2);
> end proc:

> ex:=(2*delta^2-2*beta*delta+alpha^2-2*alpha*delta+beta^2)*N^2+
  2*gamma*(beta+alpha-2*delta)*N+2*gamma^2:

> completesq(ex,[N*alpha-N*delta,N*beta-N*delta]);
                                     2                              2
           (N beta - N delta + gamma)  + (N alpha - N delta + gamma)
 

For large symbolic expressions like in this case, it may be useful showing only its structure, by replacing some (probably long) subexpressions with user-defined labels. See e.g. ?LargeExpressions,Veil .

The problem is that the product in equation (7) contains a literal subscripted name. In line print form:

proc (theta, lambda) options operator, arrow; RVS(theta, B)*(product(`&rho__l`(lambda), l = 1 .. N__op)) end proc

instead of an indexed name like in equation (4):

proc (lambda) options operator, arrow; product(rho[l](lambda), l = 6 .. N__op) end proc

This is a 2D-input trap (though there is a hint, note the difference in colour): the literal subscript is not computational as the index in rho[l] so that you cannot sum over it. This kind of trap can be avoided by using 1D-input.

@Alex Smith 

This limitation in limit arises in this line of `limit/multi` (in Maple 17):

> showstat(`limit/multi`,36..37);
`limit/multi` := proc(expr, pt, dir)
local dr, i, inds, npts, pts, r, spat, u, v, res;
       ...
  36   if not type(r,ratpoly(anything,inds)) then
  37     return FAIL
       end if;
       ...

It says that it fails for any expression that is not a rational expression (here r is the expression changed to internal variables in the set inds). Given this limitation, the problem is with ?limit,multi as it states that the type of its first argument is algebraic expression, rather than rational.

Certainly, the code in this procedure looks quite old. For example, it uses the deprecated traperror - lasterror construct. No doubt that it could be easily improved and updated. Note, on the other hand, that ?updates,Maple17,AdvancedMath#BivariateLimits describes some recent enhancements, though restricted to bivariate rational expressions.

The OP's question is about the display of a fraction, rather than the computation of a fraction which is unavoidably affected by automatic simplification rules. Now, if a pair of expressions (C,D) were handled, one for computation (as usual) and another one for display, even when C is a fixed representative within a class of mathematically equivalent expressions set by these automatic simplification rules, it could become possible to "compute" D as the display expression associated with any computational expression of this equivalence class, so that it displays in the desired form.

For pure display purposes, a demonstration procedure as in this sketch may be enough in the Standard GUI:

An additional step, copy-pasting the computational expression was described here by using the undocumented procedure Typesetting:-mcomplete, so that a variation of the above demonstration procedure could be used for the generation of its argument:

>  Frac2:=proc(x,y) sprintf("mfrac(mi(\"%A\"),mi(\"%A\"))", x,y); end proc: >  Typesetting:-mcomplete(parse(Frac2((numer,denom)(ex))),ex);   (2) >  1/3+(1/5)*x

Here the last input line was generated by copying the output and pasting as intput.

This is a hint of the things that I think could be done by decoupling computation from display when using a pair of expressions, instead of only one.

If memory becomes a constraint for a computation, it is advisable to start using the CLI (Command Line Interface) rather than a GUI, Standard in particular.

It could be useful revising the computation that has led to this long result. Perhaps some simplifications could have been made at an earlier stage. But once at this point, my suggestion is trying to put this expression in a more compact form, for better visual identification of patterns, leaving possible simplifications that could depend on properties of the Bessel functions to a second stage.

For this purpose interface labels ( see ?interface ) turn out most convenient. Note that they are available in the Classic GUI and the CLI, but not in the Standard GUI. Next, aliases ( see ?alias ) and substitutions help shortening those long Bessel function calls. And  then collecting on parameters and substituted names can make patterns visible. As an illustration, these are some steps along this goal:

> c4e := (1/2)*((2*((N*a*BesselK(1, beta)*beta*(p+pEL)*BesselI(0,
> beta*a/lambda)+N*a*BesselI(1, beta)*beta*(p+pEL)*BesselK(0,
> beta*a/lambda)+(1/2)*DD*sc*beta^2-N*a*(p+pEL))*lambda*(-1+nuf)*BesselI(1,
> beta*a/lambda)+(N*a*BesselK(1, beta)*beta*(p+pEL)*BesselI(0,
> beta*a/lambda)^2+(N*a*BesselI(1, beta)*beta*(p+pEL)*BesselK(0,
> beta*a/lambda)+(1/2)*DD*sc*beta^2-N*a*(p+pEL))*BesselI(0,
> beta*a/lambda)-(1/2)*DD*sc*BesselI(1,
> beta)*beta*(1+nuf))*a*beta))*lambda*l^2*k*BesselK(0,
> k1*lambda/l)+2*sc*N*beta^2*DD*k1*(l*lambda*BesselK(1,
> k1*lambda/l)*(-1+nuf)*BesselI(1, beta*a/lambda)+l*BesselI(0,
> beta*a/lambda)*beta*a*BesselK(1, k1*lambda/l)+BesselI(1,
> beta)*(-l*(-1+nuf)*BesselK(1, k1*a/l)+BesselK(0,
> k1*a/l)*a*k1)*lambda))/(N*beta^2*(-l*lambda*BesselI(1,
> beta*a/lambda)*BesselK(0, k2*lambda/l)*BesselK(1,
> k1*lambda/l)*k1+l*lambda*BesselI(1, beta*a/lambda)*BesselK(1,
> k2*lambda/l)*k2*BesselK(0, k1*lambda/l)+l*nuf*lambda*BesselI(1,
> beta*a/lambda)*BesselK(0, k2*lambda/l)*BesselK(1,
> k1*lambda/l)*k1-l*nuf*lambda*BesselI(1, beta*a/lambda)*BesselK(1,
> k2*lambda/l)*k2*BesselK(0, k1*lambda/l)+l*BesselI(0,
> beta*a/lambda)*beta*a*BesselK(0, k2*lambda/l)*BesselK(1,
> k1*lambda/l)*k1-l*BesselI(0, beta*a/lambda)*beta*a*BesselK(1,
> k2*lambda/l)*k2*BesselK(0, k1*lambda/l)+lambda*BesselI(1,
> beta)*k1^2*BesselK(0, k2*lambda/l)*BesselK(0, k1*a/l)*a+lambda*BesselI(1,
> beta)*k1*BesselK(0, k2*lambda/l)*l*BesselK(1, k1*a/l)-lambda*BesselI(1,
> beta)*k1*BesselK(0, k2*lambda/l)*nuf*l*BesselK(1, k1*a/l)-lambda*BesselI(1,
> beta)*k2*BesselK(0, k1*lambda/l)*l*BesselK(1, k2*a/l)+lambda*BesselI(1,
> beta)*k2*BesselK(0, k1*lambda/l)*nuf*l*BesselK(1, k2*a/l)-lambda*BesselI(1,
> beta)*k2^2*BesselK(0, k1*lambda/l)*BesselK(0, k2*a/l)*a)*DD*k):
> 
> alias(BI=BesselI,BK=BesselK):
> c:=subs(p+pEL=pp,
> BesselI(0,beta*a/lambda)=BI0,BesselI(1,beta*a/lambda)=BI1,
> BesselK(0,beta*a/lambda)=BK0,BesselK(1,beta*a/lambda)=BK1,
> c4e):
> collect(c,[pp,lambda,N,DD,beta,l,nuf,a,k,k1,k2,kk,BI1],simplify);
((((BK(1, beta) BI0 + BI(1, beta) BK0) BI1 a nuf
     + (-BK(1, beta) BI0 - BI(1, beta) BK0) BI1 a) beta - a BI1 nuf + a BI1) N
                                                       2     2    2
    lambda + (BI0 (BK(1, beta) BI0 + BI(1, beta) BK0) a  beta  - a  BI0 beta) N
              2       k1 lambda       / /      2 //
    ) lambda l  BK(0, ---------) pp  /  |N beta  ||
                          l         /   \        \\
                                                                           /
    (((%4 - %3) k1 + (%1 - %2) k2) nuf + (-%4 + %3) k1 + (%2 - %1) k2) l + |
                                                                           \
                  2       k2 lambda        k1 a
    BI(1, beta) k1  BK(0, ---------) BK(0, ----)
                              l             l
                     2       k1 lambda        k2 a \  \          /
     - BI(1, beta) k2  BK(0, ---------) BK(0, ----)| a| lambda + |
                                 l             l   /  /          \
              k2 lambda        k1 lambda
    BI0 BK(0, ---------) BK(1, ---------) k1
                  l                l
                 k2 lambda           k1 lambda \         \   \       /
     - BI0 BK(1, ---------) k2 BK(0, ---------)| a l beta| DD| + 1/2 |
                     l                   l     /         /   /       \
    /               k1 lambda                       k1 lambda \  2     2
    |sc BI1 k BK(0, ---------) nuf - sc BI1 k BK(0, ---------)| l  beta  DD
    \                   l                               l     /
          2   ///
    lambda  + |||
              \\\
    /           k1 lambda                               k1 a \
    |2 sc BK(1, ---------) BI1 - 2 sc BI(1, beta) BK(1, ----)| k1 nuf
    \               l                                    l   /
       /            k1 lambda                               k1 a \   \
     + |-2 sc BK(1, ---------) BI1 + 2 sc BI(1, beta) BK(1, ----)| k1| l
       \                l                                    l   /   /
              2                   k1 a   \     2        /
     + 2 sc k1  BI(1, beta) BK(0, ----) a| beta  DD N + |
                                   l     /              \
              2         k1 lambda      3   /
    sc BI0 a l  k BK(0, ---------) beta  + |
                            l              \
                              k1 lambda
    -sc BI(1, beta) a k BK(0, ---------) nuf
                                  l
                                k1 lambda \  2     2\   \
     - sc BI(1, beta) a k BK(0, ---------)| l  beta | DD| lambda
                                    l     /         /   /
                  3                     k1 lambda \   / /      2 //
     + 2 sc N beta  DD k1 l BI0 a BK(1, ---------)|  /  |N beta  ||
                                            l     / /   \        \\
                                                                           /
    (((%4 - %3) k1 + (%1 - %2) k2) nuf + (-%4 + %3) k1 + (%2 - %1) k2) l + |
                                                                           \
                  2       k2 lambda        k1 a
    BI(1, beta) k1  BK(0, ---------) BK(0, ----)
                              l             l
                     2       k1 lambda        k2 a \  \          /
     - BI(1, beta) k2  BK(0, ---------) BK(0, ----)| a| lambda + |
                                 l             l   /  /          \
              k2 lambda        k1 lambda
    BI0 BK(0, ---------) BK(1, ---------) k1
                  l                l
                 k2 lambda           k1 lambda \         \     \
     - BI0 BK(1, ---------) k2 BK(0, ---------)| a l beta| DD k|
                     l                   l     /         /     /
                        k1 lambda        k2 a
%1 := BI(1, beta) BK(0, ---------) BK(1, ----)
                            l             l
                k2 lambda        k1 lambda
%2 := BI1 BK(1, ---------) BK(0, ---------)
                    l                l
                        k2 lambda        k1 a
%3 := BI(1, beta) BK(0, ---------) BK(1, ----)
                            l             l
                k2 lambda        k1 lambda
%4 := BI1 BK(0, ---------) BK(1, ---------)
                    l                l

The command solve has the option DropMultiplicity for removing repeated roots. Presumably, it preserves the order, cf.:

> solve(x^2*(x-a-1)*(x-a+1),[x]);
                 [[x = 0], [x = 0], [x = 1 + a], [x = -1 + a]]
> {op(%)};
                      {[x = 0], [x = -1 + a], [x = 1 + a]}
> solve(x^2*(x-a-1)*(x-a+1),[x],DropMultiplicity=true);
                      [[x = 0], [x = 1 + a], [x = -1 + a]]

It is available since Maple 13.

@Andriy 

These two new simple problems can be handled by a simple generalization of the rule presented by Markiyan (as a conditional rule):

r:=conditional(abs(A::algebraic)^(n::integer) = 
(A*conjugate(A))^(n/2),_type(n,even)):

z:=abs(a)^2*b+abs(a)^2*c*conjugate(a);
                                     2          2   _
                           z := | a |  b + | a |  c a
applyrule(r, z);
                                   _       _2
                                 a a b + a a  c
z:=a^4*conjugate(a^4);
                                         4 _4
                                   z := a  a
z:=simplify(z);
                                            8
                                  z := | a |
applyrule(r,z);
                                      4 _4
                                     a  a
z := (a*b)^2*conjugate((a*b)^2);
                                          _______
                                     2  2   2  2
                               z := a  b  (a  b )
z := simplify(z);
                                         4  4
                                 z := | a  b  |
applyrule(r,expand(z));
                                   2  2 _2 _2
                                  a  b  a  b

Note that more complex expressions may require some additional generalization and/or additional transformation rules, depending on the algebra involved.

@Carl Love 

Debugging shows that this evalhf error appears in a procedure defined on-the-flight, specific for the system, named f and called from `DEtools/DEplot/FieldEvalAndMax2` line 9. See below the relevant excerpt of the output (where for convenience the code was run in the CLI), where the debugger is stopping at the first occurence of this error. And I am showing also the arguments to this procedure f. The third one, pt, is actually an array (an old array), while the fourth is an Array (the rtable one):

...
> stoperror(`unsupported type`):
> gradV:=map(q->simplify(diff(V,q)),[x,y]): ware float evaluation
> gradsys:=zip((a,b)->diff(a(t),t)=subs(x=x(t),y=y(t),b),[x,y],gradV):
> op(gradsys):
> display([
>    DEplot(gradsys,[x(t),y(t)],t=-20..3,
>        [seq([x(0)=2*cos(theta),y(0)=2*sin(theta)],
>            theta=[seq(k*Pi/12,k=1..23),seq(k*Pi/36,k=-2..2)])],
>        stepsize=.2,
>        title=`Gradient system`,
>        titlefont=[TIMES,BOLD,16],
>        color=blue,
>        linecolor=red,
>        arrows=MEDIUM,
>        method=rkf45),
>    cplot
>    ],view=[-2..2,-2..2]);
memory used=170.6MB, alloc=105.5MB, time=3.04
memory used=275.6MB, alloc=105.5MB, time=4.15
Error, unsupported type `%1` in evalhf
f:
   1   YP[1] := -exp(-Y[1]-1)*LambertW(exp(-1))*...
DBG> showexception
[f, "unsupported type `%1` in evalhf", complex(float)]
f:
   1   YP[1] := -exp(-Y[1]-1)*LambertW(exp(-1))*...
DBG> showstat(procname)
f := proc(N, T, Y, YP)
   1 ! YP[1] := -exp(-Y[1]-1)*LambertW(exp(-1))*(exp(-Y[1]-1)*LambertW(exp(-1))-
Y[1]*exp(-1)-LambertW(exp(-1))*exp(-1))/(exp(-1)^2+2*exp(-1)*exp(-Y[1]-1)*
LambertW(exp(-1))+exp(-Y[1]-1)^2*LambertW(exp(-1))^2);
   2   YP[2] := 1/2/(Y[2]^2+1)^(1/2)*Y[2]
end proc
DBG> args
2,
0,
pt,
Array(1..2, [0.+0.*I,0.+0.*I], datatype = complex[8])
f:
   1   YP[1] := -exp(-Y[1]-1)*LambertW(exp(-1))*...
DBG> eval(args[3])
memory used=387.9MB, alloc=121.5MB, time=5.22
array(1 .. 2,[(1)=-5.13569288283891634,(2)=-3.40084505501195711])
...

Let us look at the line 1 of f. Executed at top level it goes like:

> Y:=array(1..2,[-5.13569288283891634,-3.40084505501195711]):
> YP:=Array(1..2, [0.+0.*I,0.+0.*I], datatype = complex[8]):
> YP[1] := exp(-Y[1]-1)*LambertW(exp(-1))*...
Error, unsupported type `complex(float)` in evalhf

This error message occurs in simple assignments with an element of the Array YP in the lhs, like:

> YP[1] := LambertW(1);
Error, unsupported type `complex(float)` in evalhf

but not here when the lhs is a symbol:

 > a := LambertW(1);
                                a := LambertW(1)

So, the origin of this problem seems to lay at the "evaluation" of the lhs, when it is an element of an Array declared of datatype = complex[8], as it triggers the evaluation of the rhs using evalhf, and it causes a problem when this rhs contains LambertW, as shown before:

> Z:=Array(1..2, [0.+0.*I,0.+0.*I]):
> Z[1]:= LambertW(1);
                              Z[1] := LambertW(1)
> 
> Z:=Array(1..2, [0,0], datatype = complex[8]):
> Z[1]:= LambertW(1);
Error, unsupported type `complex(float)` in evalhf
> 
> Z:=Array(1..2, [0,0], datatype = complex[8]):
> Z[1]:= exp(1);
                                 Z[1] := exp(1)

Note that this mechanism seems independent of the value of Digits or UseHardwareFloats.

I think that this is a multifactorial regression bug. In particular, the hardware float evaluation of the rhs of the assignment without a fallback to software float evaluation seems to me a design bug. And the mix of table-arrays with rtable-Arrays at this epoch seems to me odd, if not risky. The same advice given here and elsewhere to so many users of updating their code by using rtable data structures should apply to developers also.

So, probably, workarounds should involve avoiding all the components of this bug coming together. An example is acer's suggested substitution of alpha=evalf(alpha), which avoids calling that buggy kernel routine "evalfh/LambertW".

If you want to port that MMA code into Maple by writing a procedure (or more than one), procedure local variables are probably what you need. You may initialize them in their same declaration, if you want, see ?ProgrammingGuide,Chapter06#LocalVariables . About initializing Vectors you could use their initialization syntax, something like in this toy example:

> f:=proc()
> local a:=<1,2>,b:=a[1]+1:
> b:
> end proc:
> 
> f();
                                       2

Note that use involves binding equations, but it is resolved during automatic simplification, meaning that there is no evaluation of their rhs ( cf. ?use ), so that it might not be what you need. Also, your code shows nested With statements, while nested use statements have, I think, different behavior (see examples 2 and 10 in ?use ).

You can stop at the procedure WARNING, and then look at the call stack from where it comes the call. Something like:

> stopat(WARNING):
> int(1/x, x = a[1] .. 1,method=ftoc) assuming a[1]>0;
WARNING:
   1*  print(INTERFACE_WARN(1,msg,args[2 .. -1]))
DBG> where 2
IntegrationTools:-Definite:-DirectFinishInt: WARNING(op(t))
	["unable to determine if %1 is between %2 and %3; try to use assumptions or use the AllSolutions option", 0, a[1], 1]
WARNING:
   1*  print(INTERFACE_WARN(1,msg,args[2 .. -1]))
DBG> quit