@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.

Why not use LinearAlgebra:-Eigenvector?

Choosing y ' (0) = 0 as suggested by tomleslie works fine numerically on the interval 0..1:

restart;
edo := -b*y(x)*(diff(y(x), x))^2*x-a*y(x)*x+(diff(y(x), x, x))*x+2*(diff(y(x), x)) = 0
edo1 := subs(a = 1, b = 1, edo);
res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..1);
plots:-odeplot(res,[x,y(x)],0..1);

## To push it further, you have to work harder:

res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..2,approxsoln=[y(x)=1-0.3*x^2],initmesh=256);
plots:-odeplot(res,[x,y(x)],0..2);

The solution to the singular initial value problem {edo1, y(0)=1, D(y)(0)=0} probably tends to infinity at some finite x0, i.e. y(x)-> infinity for x-> x0 from the left. x0 seems to be greater than 2 and probably less than 3.

Note. There is some evidence for an analytic solution at zero having radius of convergence less than 2.141580437. Powers all even. Using the terms up to degree 50 of that solution (sol) as approxsoln in dsolve you can get success with
res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..2.09,approxsoln=[sol],initmesh=1024);
plots:-odeplot(res);

@ecterrab We agree that a syntax like f(a)(b) certainly is useful as in e.g. D(f)(b).
I was deliberately vague when calling f(t)(t) a weird construction.

But notice that it was the appearance of several "weird" constructions like that that led me to find Torleif's mistake, which he admits to having made, when making last minute changes to a worksheet.

So far nobody asked you why you need this function. What do you intend to use it for? The answer to that may determine which method is preferable.

@ecterrab Sure, but can you give a non-weird example of f(t)(t)? Point being that there are 2 t's.

D(sin)(sin) would not make much sense, would it?

@torabi What sense does it have to remove nonlinear terms?