EE does not contain a suboperand sqrt(Bx1^2+By1^2), thus it won't be seen by eval.

Try this:

EE:=1/sqrt(a^2+b^2)*(c+d);
op(EE);
op(%[1]);
eval(EE,1/sqrt(a^2+b^2)=1/d1);
#or instead:
eval(EE,a^2+b^2=d1^2);
simplify(%) assuming d1>0;

When using Digits=d d significant figures are carried in intermediate calculations. In Maple there is no claim that this results in d significant digits in the result, which it often does not.
To explain the behavior you mention take a look at
showstat(`evalf/frac`);

Markiyan Hirnyk already took care of the color part. As for the legend you can do as follows:

restart;
R[123]:=(1+sin(lambda*S))/2;
CL:=[red,blue,green,yellow,orange,"SkyBlue","PeachPuff"]:
for C to 7 do
S := .7+.1*C;
p || C := plot(R[123], lambda = 5 .. 25, y = 0 .. 1, legend = typeset("S=",0.7+0.1*C),color=CL[C])
end do :

plots:-display(p1,p2,p3,p4,p5,p6,p7);

(1) In the beginning you assign IT to a matrix whose elements contain names Ixx, Ixy etc.
These names appear in dsys because they were never assigned values. Much later you redefine IT to a matrix with numerical elements, but that doesn't give Ixx, Ixy, ... numeric values. One simple way out of this is to assign directly to Ixx, Ixy, etc. There are other ways of course.
(2) In your ini phiy = .. and phiz = .. are missing (0), i.e. they should be phiy(0) = .. and phiz(0) = ..
(3) In L you are using two different phiz's: try indets(L,function); You will see that phiz(t) appears twice!. Thus even after the substitutions phiz(t) is not gone from Lb. 1D input is much to be preferred in my opinion, then this couldn't have happened. Even MaplePrimes tries to prevent me from telling you how these two differ. I tried twice. Now the 3rd time: The one is  varphi;z(t)  the other  phi;z(t) , deliberately missing part of the characters this 3rd time.
(4) You assign R:=0 early on, yet later in Eq18 you substitute into R??

Here follows a linear algebra version.

restart;
expand((d^2-a^2)^2);
ode1:=diff(H(z),z$4)-2*a^2*diff(H(z),z,z)+a^4*H(z)=a^2*lambda*alpha(z);
ode2:=diff(alpha(z),z$2)-a^2*alpha(z)=-lambda*H(z);
bcs:=alpha(0)=0,alpha(1)=0,H(0)=0,H(1)=0,(D@@2)(H)(0)=0,(D@@2)(H)(1)=0;
#Set u=(D^2-a^2)H (not Maple notation). Introduce the vector function q:
q:=<H,u,alpha>;
#With L defined as
L:=Matrix([[0,1,0],[0,0,a^2*lambda],[-lambda,0,0]]);
#the system can be written (D^2-a^2)q = L.q (still not in Maple notation)
#We decouple the system:
ev,V:=LinearAlgebra:-Eigenvectors(L);
ev:=simplify(ev) assuming lambda>0,a>0;
#Find solutions for the decoupled system with p=V^(-1).q. Notice that it follows from bcs that the q's as well as the p's satisfy zero boundary conditions at 0 and 1.
#The third eigenvalue of L is real, so we consider that.
eq1:=diff(p3(z),z,z)-a^2*p3(z)=ev[3]*p3(z);
dsolve({eq1,p3(0)=0,p3(1)=0},{p3(z),lambda});
res:=simplify(%[3]) assuming positive;
odetest({H(z)=sin(n*Pi*z),alpha(z)=c*sin(n*Pi*z)},{ode1,ode2});
remove~(has,%,sin);
eval(%,lambda = (a^2+Pi^2*n^2)^(3/2)/a);
solve(%,c);

Edited: I made a simple mistake below: The coefficient of H(z) in ode1 should be a^4, not a^2.
I inserted a correction and commented out my original ode1. Funny that I should also make a mistake in the ode that I corrected in your version. This also means that the eigenvalue is wrong. It has also been corrected.
You need to correct ode1 to correspond to the image. Your ode1 as defined in the code is nonlinear but only of second order.
But the problem is an eigenvalue problem for a system of linear differential equations. You need to determine lambda (eigenvalue) so that the problem has nontrivial solutions and also determine those (eigenfunctions).

restart;
#To get the proper ode1 the easiest seems to be to expand as follows:
expand((d^2-a^2)^2);
#Then write the ode1 accordingly:
#ode1:=diff(H(z),z$4)-2*a^2*diff(H(z),z,z)+a^2*H(z)=a^2*lambda*alpha(z);
ode1:=diff(H(z),z$4)-2*a^2*diff(H(z),z,z)+a^2*H(z)=a^2*lambda*alpha(z);
ode2:=diff(alpha(z),z$2)-a^2*alpha(z)=-lambda*H(z);

bcs:=alpha(0)=0,alpha(1)=0,H(0)=0,H(1)=0,(D@@2)(H)(0)=0,(D@@2)(H)(1)=0;
#In general there is only the trivial solution:
dsolve({ode1,ode2,bcs});
#Testing your proposed solution you see the problem:
odetest({H(z)=sin(Pi*z),alpha(z)=sqrt(Pi^2+a^2)/a*sin(Pi*z)},[ode1,ode2,bcs]);
#This reveals that all 6 bcs are satisfied, but that the ode1 and ode2 are satisfied iff
lambda=(a^2+Pi^2)^(3/2)/a; #Corrected as explained above!
#Maple doesn't refuse to deal with the problem of finding also lambda. The syntax would be:
dsolve({ode1,ode2,bcs},{alpha(z),H(z),lambda});
#Actually you get a result, but it is not pretty as it stands.

Maple only does the separation of variables:
restart;
PDE:=-diff(f(x1,x2),x1,x1)-diff(f(x1,x2),x2,x2) = a^2*f(x1,x2);
pdsolve(PDE);

IBC should be

IBC := {D[1](f)(0,x2) = 0, D[1](f)(L1,x2) = 0, D[2](f)(x1,0) = 0, D[2](f)(x1,L2) = 0};

Maple can solve the ODEs:

res:=pdsolve(PDE);
odes:=op([2,1],res);
dsolve(odes[1]) assuming _c[1]<0; #if that is the relevant assumption
dsolve(odes[2]);

Your system is not linear!
But you could try

dsolve(sys,vars,type=series,order=4);

(The default order is 6).

The purely symbolic result is wrong as can already be seen from just taking the first term.
The assumptions don't help any. Initially I use Int instead of int since I want to do numerical integration as test at the end.

restart;
#assume(0 < ct, ct < cl, 0 < H, 0 < tau, 0 < omega);
u := Int(Int(sin(omega*(ts+ys/ct))*Heaviside(t-ts-ys/ct), ys = H-cl*ts .. ct*ts), ts = H/(cl+ct) .. (H+cl*tau)/(cl+ct));
#value turns Int into int:
res1:=eval(value(u),[cl=6300,ct=3140]):
res2:=value(eval(u,[cl=6300,ct=3140])):
simplify(res1-res2):
evalf(eval(%,{omega=1,t=1,tau=1,H=1}));
#We see that we don't get anything close to zero.
#This will indicate that the result res2, where the substitution is done first, is the correct one:
evalf(eval(res2,{omega=1,t=1,tau=1,H=1}));
#Here the integration is purely numerical (because of Int instead of int (and no 'value'):
evalf(eval(u,{cl=6300,ct=3140,omega=1,t=1,tau=1,H=1}));

You will have to write a procedure, but it should be fairly easy as convertsys provides the information. Here I use the first example in the help page for convertsys. I have not written a procedure, but the idea can be used.

restart;
deq1 := diff(y(t), t$2) = y(t)-x(t), diff(x(t), t) = x(t):
init1 := y(0) = 1, D(y)(0) = 2, x(0) = 3:
#It is important that you provide your own names for the output dependent variables (I used Y and YP) since otherwise they will be local variables not immediately accessible:
res:=DEtools[convertsys]({deq1}, {init1}, {x(t), y(t)}, t,Y,YP);
n:=nops(res[1]);
res12:=eval(res[1..2],{seq(Y[i]=Y[i](t),i=1..n),seq(YP[i]=diff(Y[i](t),t),i=1..n)});
eqs:=convert(res12[1],set);
ics:={seq(Y[i](res[3])=res[4][i],i=1..n)};
dsolve(eqs union ics);
eval(%,res12[2]);
remove(has,%,diff);
#Compare with the direct approach:
dsolve({deq1,init1});
#Without initial conditions:
res:=DEtools[convertsys]({deq1}, {}, {x(t), y(t)}, t,Y,YP);
n:=nops(res[1]);
res12:=eval(res[1..2],{seq(Y[i]=Y[i](t),i=1..n),seq(YP[i]=diff(Y[i](t),t),i=1..n)});
eqs:=convert(res12[1],set);
if not res[3]='undefined' then ics:={seq(Y[i](res[3])=res[4][i],i=1..n)} else ics:={} end if;
dsolve(eqs union ics);
eval(%,res12[2]);
remove(has,%,diff);
#Comparison:
dsolve({deq1});

Reading the help page for MatrixFunction (at this time in Maple 12) shows that the definition used is the one presented in Gantmacher, Matrix Theory (vol. I). The definition is given on p. 96. By only reading pp. 95-97 we see that ln(1+J) can be found like this:

r:=unapply(add(a[k]*lambda^k,k=0..2),lambda);
f:=x->ln(1+x);
eqs:={r(2)=f(2),D(r)(2)=D(f)(2),(D@D)(r)(2)=(D@D)(f)(2)}; #Gantmacher (6) p.97
solve(eqs);
R:=eval(r(lambda),%);
eval(R,lambda=J);

You could try a fixed stepsize, which entails using a classical method, e.g. rk4:

S10 := dsolve({ics, de}, numeric, method = classical[rk4],known=F,' output' = Array([seq(20/200*i, i = 0 .. 200)]));
gr:=plots:-odeplot(S10, [t, x(t)],0..20,axes=boxed,color=red,linestyle=solid,legend = ["Random force "]):gr;

However, I'm not sure I understand what you want to do.

ineq:=diff(N(t),t)-mu*N(t)<=Pi;
map(curry(`*`,exp(-mu*t)),ineq);
#For T>=0 this is OK:
map(Int,%,t=0..T);
value(%);
#Unfortunately, int doesn't immediately see what we see. However:
IntegrationTools:-Parts(%,1);
solve(%,{N(T)}) assuming real;

Since the uploaded worksheet doesn't contain the definition of Platform[i] I cannot be totally sure what the problem is, but since the animation parameters is given to Platform[i] as a float and since i in the sequence of polygonplots Platform[i] , i = 2..200, refers to integers you may get an error like

Error, (in plots/animate) expecting plot structure but received: Platform[2.]

and remember also to give the number of frames (199) since otherwise you will get the default. The remedy: Use round.

animate(display,[Platform[round(ii)]],ii=2..200,frames=199,background=B);

The use of 'ii' instead of 'i' is only necessary if 'i' is assigned, which it will be if Platform[i] has been defined in a do loop.

Another (and simpler) way is to use display/insequence:

display(seq(display(B,Platform[i]),i=2..200),insequence=true);

It is easier to help if you type the expression instead of giving us an image. Most of us would like to try if our suggestions actually work and that means typing it ourselves if you give us only a picture.

Now replace the inert Diff by diff (thanks to Axel Vogt): 
sec1:=-2*I*sqrt(beta)*diff(p(T)+I*q(T),T)+(1-beta)*(p(T)-I*q(T))*(r(T)+I*s(T))*exp(I*sigma*T)-delta*I*sqrt(beta)*(p(T)+I*q(T))=0;
evalc(sec1) assuming beta>0;