Could you clarify?

while 1=2 do blah end do;

Change Epi to:

Epi := proc (t) display(circle([R*cos(omega*t), R*sin(omega*t)], r, color = red), disk([R*cos(omega*t)+r*cos(k*t), R*sin(omega*t)+r*sin(k*t)], .2, color = blue),plot([R*cos(omega*s)+r*cos(k*s), R*sin(omega*s)+r*sin(k*s),s=0..t])) end proc;

I tried saving the worksheet with output and using 100 frames. On my 5 year old computer it took just a second or two. (Otherwise I never save worksheets with output!)

If the matrix or vector is called A, then you can convert it to a list ( convert(A,list) ), then sort that using the sorting procedure (s,t)->evalb(abs(s)<abs(t))

sort(convert(A,list),(s,t)->evalb(abs(s)<abs(t)));

#If A is actually a vector then you don't need to convert to a list:

sort(A,(s,t)->evalb(abs(s)<abs(t)));

The following seems to me to suggest that the integral may be exactly 0.
For convenience the expression has been slightly rewritten.

f221 := 0.1255914626e-2*375*(sin(2.*a)*cos(b)*cos(w)-0.4)*Heaviside(abs(sin(2.*a)*cos(b)*cos(w)-sin(a)^2*sin(b)*sin(w))-0.4)*sin(a)*cos(b)^3*sin(w)*(1.+.25*cos(2.*b));
with(plots):
animate(plot3d,[f221,a=0..2*Pi,w=0..2*Pi,axes=boxed],b=0..Pi/2);
h1,f:=selectremove(type,f221,specfunc(anything,Heaviside));
h:=op(h1);
animate(plot3d,[f,a=0..2*Pi,w=0..2*Pi,axes=boxed],b=0..Pi/2);
#Without the Heaviside factor the integral is 0:
int(f,[a=0..2*Pi,w=0..2*Pi,b=0..Pi/2]);
#The symmetry exhibited here suggests that possibly the integral is zero:
animate(implicitplot,[h=0,a=0..2*Pi,w=0..2*Pi,gridrefine=1],b=0..Pi/2);

Added: Replacing Heaviside with an approximation enforces the suspicion that the result is exactly zero:

H:=(c,x)->1/2+tanh(c*x)/2;
infolevel[`evalf/int`]:=4:
evalf[15](Int(f*H(5,h),[a=0..2*Pi,w=0..2*Pi,b=0..Pi/2],epsilon=0.1));
evalf[15](Int(f*H(10,h),[a=0..2*Pi,w=0..2*Pi,b=0..Pi/2],epsilon=0.1));
evalf[15](Int(f*H(20,h),[a=0..2*Pi,w=0..2*Pi,b=0..Pi/2],epsilon=0.1));
evalf[15](Int(f*H(40,h),[a=0..2*Pi,w=0..2*Pi,b=0..Pi/2],epsilon=0.1));
evalf[15](Int(f*H(50,h),[a=0..2*Pi,w=0..2*Pi,b=0..Pi/2],epsilon=0.1));

Does the following do what you want?

param:=convert(indets(sols,name),list);
M:=Matrix([param,seq(subs(k,param),k=sols)]);
ExportMatrix("F:/MapleDiverse/MaplePrimes12-10-11.txt",M);
restart;
ImportMatrix("F:/MapleDiverse/MaplePrimes12-10-11.txt");

You may also try solve:

res:=solve(sol1=0,t);

#with result

res := (50/49)*ln(exp(56/5)-sqrt((exp(56/5))^2-1))*sqrt(35), (50/49)*ln(exp(56/5)+sqrt((exp(56/5))^2-1))*sqrt(35)

evalf[20](res);

Notice the use of a higher value for Digits due to the huge and nearly equal terms in the logarithm of the first solution.

The second question: Time to obtain velocity v0:

solve(diff(sol1,t)=v0,t);

I hope I got it right. I removed all that have to do with the deprecated linalg package and also removed unnecessary conversions. It doesn't matter in this case, but the constant pi is spelled Pi in Maple (capital P).
The result is 8 and comes almost instantaneously.

restart;
#interface(rtablesize=12):
a1:=1:
a2:=1:
T:=[x1,x2,x3,x4,y1,y2,y3,y4,u1,u2,u3,u4]: #List
dx:=sqrt((x1-x3)^2+(x2-x4)^2):
x1dot:=a2/(2*Pi*dx^2)*<x4-x2,x1-x3>:
#Pi not pi
x2dot:=a1/(2*Pi*dx^2)*<x2-x4,x3-x1>:
dy1:=sqrt((x1-y1)^2+(x2-y2)^2):
dy2:=sqrt((x3-y1)^2+(x4-y2)^2):
dy3:=sqrt((x3-y3)^2+(x4-y4)^2):
dy4:=sqrt((x3-y3)^2+(x4-y4)^2):
y1dot:=a1/(2*Pi*dy1^2)*<x2-y2,y1-x1>+a2/(2*Pi*dy2^2)*<x4-y2,y1-x3>+<u1,u2>:
y2dot:=a1/(2*Pi*dy3^2)*<x2-y4,y3-x1>+a2/(2*Pi*dy4^2)*<x4-y4,y3-x3>+<u3,u4>:
xdot:=<x1dot,x2dot>:
zdot:=<xdot,y1dot,y2dot>:
l:=<zdot,-u1,u2,-u3,u4>:
g1:=<0,0,0,0,0,0,0,0,1,0,0,0>:
g2:=<0,0,0,0,0,0,0,0,0,1,0,0>:
g3:=<0,0,0,0,0,0,0,0,0,0,1,0>:
g4:=<0,0,0,0,0,0,0,0,0,0,0,1>:
#Lie Bracket Function
LB:=proc(m::Vector,w::Vector)
 VectorCalculus:-Jacobian(w,T).m-VectorCalculus:-Jacobian(m,T).w;
end proc;
C:=Matrix([g1,g2,g3,g4,LB(l,g1),LB(l,g2),LB(l,g3),LB(l,g4),LB(l,LB(l,g1)),LB(l,LB(l,g2)),LB(l,LB(l,g3)),LB(l,LB(l,g4))]):
LinearAlgebra:-Rank(C);

I wonder why you would want to do it, but certainly getdata will get you the data:

contourplot(x*F1,x=-3..3,eta=-0.2..0.2);
dat:=plottools:-getdata(%):
dm:=op(op~(3,[dat])):
M:=<dm>;
#Notice that if you plot M by plot(M); you get connecting lines unless you ask for style=point, whereas plotting [dm] doesn't.
plot([dm]);
#You can use ExportMatrix to export M and ImportMatrix to get it back. The sequence dm is easily recovered from M:
dm[1..10];
seq(M[i..i+1,1..2],i=1..20,2);

restart;
_EnvFormal:=false;
sum(log(n)/n,n=1..infinity);
?sum,details

Added: You talk about the relevance. Well, the error can be reproduced like this:

restart;
sum(log(n)/n^p,n=1..infinity);
eval(%,p=1);

Once you notice that implicitplot doesn't find anything outside a certain rectangle, you can redraw using that rectangle. You could try (first perhaps without the optional time consuming additions):

animate(implicitplot,[p(N, 0.5, s, 0.1,-1, f2) =0, s=-5..5, f2 = -5..5, color = black,grid = [35, 35], gridrefine = 2, resolution = 1000, crossingrefine = 3],N=1.35..3,frames=10);

You will see something interesting taking place between N = 1.5 and N=1.9.

Thus you could continue with

animate(implicitplot,[p(N, 0.5, s, 0.1,-1, f2) =0, s=-5..5, f2 = -5..5, color = black,grid = [35, 35], gridrefine = 2, resolution = 1000, crossingrefine = 3],N=1.5..1.9,frames=10);

or just use more frames.

Edited: Did you mean to arrange the parameters in the order you do:

In dsolve: [M, s, k, f1,f2] and in p: proc(N, s, M, k,f1,f2) ?

I tried the code and had absolutely no problem. Have you tried after a restart? Or have you omitted some lines?

frequency := [1 = 6, 2 = 7, 3 = 91, 4 = 10];
(`[]`@op)~(frequency);

#To understand what is going on you can take a look at the outputs of the following two commands. f~ means elementwise use of f and @ is used for composition.
op~(frequency);
`[]`~(frequency);

A version more easily understood perhaps is

map(eq->[lhs(eq),rhs(eq)], frequency);

_Z is just the name used in RootOf for the unknown. The answer from solve in this case is really just a rewriting of the question. This may seem silly, but Maple can work with RootOf expressions in several ways.

Besides that, you should use exp(-x) instead of e^(-x) if you are using the exponential function.

I have been sorting the output according to the sizes of the eigenvalues using the following procedure, which also contains a description:

EigenSort:=proc(ev::Vector,P::Matrix,
krit::procedure:=((s,t)->evalb(
evalf(Re(s)<Re(t)) or evalf(Re(s)=Re(t)) and evalf(Im(s)>=Im(t)) )))
local n,EV,SEV,L,krit2;
description "Sorts output from Eigenvectors(A). As an optional third argument can be given the sorting procedure. The default procedure sorts according to the real parts of the eigenvalues. For a pair of complex conjugates the one with positive imaginary part is taken first.
The default procedure uses evalf before compairing.";
uses LinearAlgebra;
n:=Dimension(ev);
EV:=convert(<ev|<seq(i,i=1..n)>>,listlist);
krit2:=(s,t)->krit(op(1,s),op(1,t));
try
SEV:=sort(EV,krit2);
catch:
error cat(procname," cannot sort the given input");
end try;
L:=map2(op,2,SEV);
<seq(ev[L[i]],i=1..n)>,Matrix([seq(Column(P,L[i]),i=1..n)]);
end proc:


Describe(EigenSort);
A:=LinearAlgebra:-RandomMatrix(2);
EV1:=LinearAlgebra:-Eigenvectors(A);
EV2:=EigenSort(EV1);

Assuming that the image part of your code is

N:=t->( (1/N0-1)*exp(-2*r*t)+1)^(-1/2);

then N0 =1 would give you the constant function 1. Also try other values of N0.