What I see is (for the initial value problem, as an example):

y(x) = JacobiSN((1/2)*sqrt(2)*x+InverseJacobiSN(1, I), I)

and after a % I get

y(x) = JacobiSN((1/2)*sqrt(2)*x+EllipticK(I), I)

How did you get those sn's in the first place? Is there in fact something missing in your worksheet?

@Markiyan Hirnyk Since you need a restart in between, it is difficult to compare two exact result that are so big:

But just try

restart;
int(......);
evalf(%);

restart;
int(......);
evalf(%);

There you surely don't have numerical integration, so the cause is roundoff.

@Konstantin@ Yes, if int cannot do the integration (as is the case here), then it uses numerical integration.

Try

int(exp(sin(x)+x+exp(x)),x=0..1);

@Bendesarts The code I used:

restart;
K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);
omega[sw]:=beta/(1-beta)*omega[s];
for i to 4
do
r[i]:=sqrt((u[i](t)/(L/2))^2+(v[i](t)/H)^2):
omega[i]:=omega[st]/(1+exp(b*v[i](t)))+omega[sw]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)+omega[i]*(L/2)/H*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*(L/2)/H*v[i](t)+(K.<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do;
paramsGeo:=L=0.015,H=0.015,beta=0.5,Vf=0.3;
omegaS:=eval(Pi*Vf/L, [paramsGeo]);
paramsCycle:=omega[s]=omegaS,Au=1,Av=1,b=100;
params:=paramsGeo,paramsCycle;

sys:=map(op,eval({seq(EqSys[i],i=1..4)},{params}));
#ic:={u[1](0)=0.8, v[1](0)=0,u[2](0)=0.8, v[2](0)=0,u[3](0)=0.8, v[3](0)=0,u[4](0)=0.8, v[4](0)=0};
ic:={u[1](0)=0.8, v[1](0)=0.1,u[2](0)=0.8, v[2](0)=0.1,u[3](0)=0.8, v[3](0)=0.1,u[4](0)=0.8, v[4](0)=0.1};
res:=dsolve(eval(sys,omega[st]=50) union ic,numeric);
plots:-odeplot(res,[u[1](t),v[1](t)],0..1);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..1,thickness=2);
#The last plot:

I don't see any text in the body of the question. How come Carl can see it?

@GambiaMan What is it you want plotted: Only the points corresponding to the times in your sample? If that is the case your option numpoints=2000 confuses me.

@Bendesarts Change v to

vt:=<seq(a[i](t),i=1..4)>;

Without changing v:

v:=Vector(4,symbol=a);
vt:=apply~(v,t);

Comment: Why not just v(t)? Why this funny apply, which is otherwise almost never used?
Because () as well as [] are used these days for getting elements from vectors and matrices.
Thus v(1) returns a[1], and v(t) is understood as the t'th element of v, which doesn't make sense if t is not an integer.
There are differences between () and [], though. Look for programmer indexing contra mathematical indexing.
See
?Indexing Arrays, Matrices, and Vectors

## Making a function V which returns vectors:
V:=t->apply~(v,t);
V(s);
##Or the slightly strange looking version of the latter:
V:=unapply(apply~(v,t),t);
V(6);

@red1eco With the values p=50, q=0.01, W=0.1, I'm not surprised to find no visible change in the graph of x when tau is varied from 0 to 0.2.
In fact I made an animation in the parameter tau from 0 to 0.2. This animation wasn't very animated. For the fun of it I tried tau=0..20. That was much more exciting.
##At the bottom I have looked at differences instead. For that purpose I added the option output=listprocedure in dsolve below.
Note: In the first version I used the name p for the procedure I now call P (capital P). Since p is one of the parameters that had to be changed.

In your last worksheet to show the different graphs in the same picture you just need to save each plot in a variable right after it is produced, like this (e.g.) for the first one

plots:-odeplot(res, [t, x(t)], 0 .. 9, linestyle = dot);  p1 := %:

If you then have made the plots p1,p2,p3,p4,p5 you can display them together using display from the plots package:

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

You will be disappointed, however, as I warned you above: they fall almost on top of each other.

##########
The animation and the parameters option:
restart;
sys2 := [diff(x(t), t) = r*x(t)*(K-x(t)-a*y(t))/K, diff(y(t), t) = s*y(t)*(L-y(t)+b*x(t))/L-q*y(t)*E(t), diff(E(t), t) = W*((p-tau)*q*y(t)-c)*E(t)];
params := convert(indets(sys2, name) minus {t}, list); #LIST
paramvalsFix := {K = 2000, L = 1200, W = .1, a = 1.2, b = .2, c = 100, p = 50, q = 0.1e-1, r = 2, s = 1.2};
#tau left out above since we want to vary that
res:=dsolve({sys2[], E(0) = 75, x(0) = 451, y(0) = 1290}, numeric,parameters=params,output=listprocedure);
res(parameters=[op(paramvalsFix),tau=0]); #Setting parameter values
plots:-odeplot(res, [t, x(t)], 0 .. 9, linestyle = dot); #Example plot
## Now setting up a procedure to do the parameter setting and to produce the plot.
##The procedure allows optional arguments and those will be passed to odeplot (the _rest inside P).
P:=proc(tt) res(parameters=[op(paramvalsFix),tau=tt]); plots:-odeplot(res, [t, x(t)], 0 .. 9,_rest) end proc;
P(0); #Test with tau=0
P(0,color=blue); #Test with tau=0 and using the color blue
##Animation using p:
plots:-animate(P,[tau,caption="Superply dull"],tau=0..0.2);
plots:-animate(P,[tau,caption="More exciting and blue",color=blue],tau=0..20);
##The following indicates how you can display several graphs of x using the procedure P:
plots:-display(P(0,linestyle=dot),P(5,linestyle=dash),P(20,color=blue));
########DIFFERENCES:
X:=subs(res,x(t)); #Saving the procedure for computing x(t) in X (capital X).
N:=200: #Number of t-points
T:=<seq(9.0*i/N,i=0..N)>; #t-point vector
T[-1]; #The last point
##Now a procedure that produces a corresponding vector of x-values:
Q:=proc(tt) res(parameters=[op(paramvalsFix),tau=tt]); X~(T) end proc;
plot(T,Q(0)-Q(0.05),caption="The difference between x values for tau=0 and tau=0.05");

Maybe replacing sqrt(Drag*g/m) by sqrt(g/m)*drag2. The first appearance of Drag should then be replaced by drag2^2.

I have no reasonable way of testing it, since I don't have your data.

Imaginary values could appear this way:

ln(cosh(sqrt(-65.)));
                        -1.576130914+3.141592654*I

@red1eco In order to continue you can then do

paramvals := {K = 2000, L = 1200, W = .1, a = 1.2, b = .2, c = 100, p = 50, q = 0.1e-1, r = 2, s = 1.2, tau = .1};
res:=dsolve(eval({sys2[],x(0)=1.2,y(0)=0.5,E(0)=.1},paramvals),numeric);
plots:-odeplot(res,[t,x(t)],0..1);
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,50*E(t)]],0..1,thickness=3,legend=[x,y,E]);


The same initial values have been used.

(I use my real name as user name. I retired from the Technical University of Denmark in 2011.)

@acer To compare the 3 methods I tried with an ode having a symbolic solution:

ode:=diff(x(t), t, t)+diff(x(t), t)+x(t) = 0;

Using the same initial values x(0)=-1,D(x)(0)=-0.8 and evaluating the second derivative at t = 4*Pi/sqrt(3).

The exact solution is easily found by

res:=dsolve({ode,x(0)=-1,D(x)(0)=-0.8});

#to be res:=x(t) = -13/15*sqrt(3)*exp(-1/2*t)*sin(1/2*sqrt(3)*t)-exp(-1/2*t)*cos(1/2*sqrt(3)*t)
##The result for x '' (t) at  t = 4*Pi/sqrt(3) is
rE:=eval(diff(rhs(res),t,t),t=4*Pi/sqrt(3));
r4:=evalf(rE):

The differences r1-r4, r2-r4, r3-r4 for the 3 methods (with Digits=10) were:

-0.301619e-5, -8.60919101217106*10^(-8), 9.21698134481730*10^(-12)

I also tried Digits = 15. That was disastrous for the plot in method 1 and the results were:
-.289843880257555, -8.61594650208852*10^(-8), -5.83379178298316*10^(-11)

It has been pointed out (coming indirectly from Allan Wittkopf, I think) that in general it is best to let dsolve/numeric do as much of the computation as possible.
This certainly confirms that.

@necron When I download from the link I still get t1.
However, I changed t1 to t myself and had no luck: In fact the system lost contact with the kernel when I tried your worksheet.

You are trying to evaluate the double integral numerically, but what is t1?

@red1eco You need concrete values for the parameters. I chose them randomly as positive floating point numbers between 1 and 9. You also need initial values. I took them out of my head:
restart;
sys2 := [diff(x(t), t) = r*x(t)*(K-x(t)-a*y(t))/K, diff(y(t), t) = s*y(t)*(L-y(t)+b*x(t))/L-q*y(t)*E(t), diff(E(t), t) = W*((p-tau)*l*y(t)-c)*E(t)];
params := `minus`(indets(sys2, name), {t});
r:=rand(1..9.0);
r(); #Test of the random number generator r
paramsR:={seq(p=r(),p=params)};
#Your sys2 is a list. By sys2[] you get the sequence of elements in the list:
res:=dsolve(eval({sys2[],x(0)=1.2,y(0)=0.5,E(0)=.1},paramsR),numeric);
plots:-odeplot(res,[t,x(t)],0..1); #Just the graph of x
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,E(t)]],0..1,thickness=3,legend=[x,y,E]); #All in the same plot

@Markiyan Hirnyk Markiyan, I'm not about to enter that old discussion with you about stability of equilibria for nonlinear autonomous systems. These results are extremely well known by people who have had a decent course in differential equations. Since I know that you are Ukrainian, I believe that you speak Russian (I don't). There are quite a few very good textbooks on odes written by Russians. I know that because they have also been translated into English.