You must have meant

E, V:= Eigenvectors(A);

Information about the known option is given in the help page

?dsolve,numeric

Here is what you can do in the example you give.

restart;
u := 8.53;
dsys0:={diff(q(t), t,t) = u*(1 - q(t)^2) * diff(q(t),t) - q(t), q(0) = 0, D(q)(0) = 1};
dsol0:=dsolve(dsys0,numeric,output=listprocedure);

yy0,yy1:=op(eval([q(t),diff(q(t),t)],dsol0));
#In your equation dsys10 you replace q(t) by yy0(t) and diff(q(t),t) by yy1(t) as I have done here.
dsys10:={q10(0) = 0, diff(q10(t), t, t) = (8.53*(1-yy0(t)^2))*(diff(q10(t), t))-17.06*yy0(t)*yy1(t)*(diff(q10(t), t))-q10(t), (D(q10))(0) = 1};
dsol10:=dsolve(dsys10,numeric,output=listprocedure,known=[yy0,yy1]);
#Plotting just to see that it works.
plots:-odeplot(dsol10,[t,q10(t)],0..5);

Information about the known option is given in the help page

?dsolve,numeric

Here is what you can do in the example you give.

restart;
u := 8.53;
dsys0:={diff(q(t), t,t) = u*(1 - q(t)^2) * diff(q(t),t) - q(t), q(0) = 0, D(q)(0) = 1};
dsol0:=dsolve(dsys0,numeric,output=listprocedure);

yy0,yy1:=op(eval([q(t),diff(q(t),t)],dsol0));
#In your equation dsys10 you replace q(t) by yy0(t) and diff(q(t),t) by yy1(t) as I have done here.
dsys10:={q10(0) = 0, diff(q10(t), t, t) = (8.53*(1-yy0(t)^2))*(diff(q10(t), t))-17.06*yy0(t)*yy1(t)*(diff(q10(t), t))-q10(t), (D(q10))(0) = 1};
dsol10:=dsolve(dsys10,numeric,output=listprocedure,known=[yy0,yy1]);
#Plotting just to see that it works.
plots:-odeplot(dsol10,[t,q10(t)],0..5);

@goli The problem is that the starting guess for z = 0 is H = 0 and you are dividing by H^2 in eq(z).

You can modify the starting guess to e.g. H=max(h*sqrt(z^3*m),h) as in the following.

restart;
eq := z-> 1=(m*(1+z)^3-k*(1+z)^2)*h^2/(H^2)+(((2*n-1)-(k*(1+z)^2*h^2/H^2))*((1-m+k)/(2*n-1+k))*((((H^2/h^2)-k*(1+z)^2)/(1-k))^(n-1)));

Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=max(h*sqrt(z^3*m),h) ) end if;
m := 0.34: k := 0.04: n := -0.34: h := 0.69:
l := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)*h/Y(z),L(0)=0}, type=numeric,range=0..2, known=Y);
A1 := ((0.35/Y(0.35))*(int(1/Y(z), z = 0 .. 0.35))^2)^(1/3);   
R := h*evalf(sqrt(m)*(int(1/Y(z), z = 0 .. 1091.3)));

@goli The problem is that the starting guess for z = 0 is H = 0 and you are dividing by H^2 in eq(z).

You can modify the starting guess to e.g. H=max(h*sqrt(z^3*m),h) as in the following.

restart;
eq := z-> 1=(m*(1+z)^3-k*(1+z)^2)*h^2/(H^2)+(((2*n-1)-(k*(1+z)^2*h^2/H^2))*((1-m+k)/(2*n-1+k))*((((H^2/h^2)-k*(1+z)^2)/(1-k))^(n-1)));

Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=max(h*sqrt(z^3*m),h) ) end if;
m := 0.34: k := 0.04: n := -0.34: h := 0.69:
l := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)*h/Y(z),L(0)=0}, type=numeric,range=0..2, known=Y);
A1 := ((0.35/Y(0.35))*(int(1/Y(z), z = 0 .. 0.35))^2)^(1/3);   
R := h*evalf(sqrt(m)*(int(1/Y(z), z = 0 .. 1091.3)));

@goli For large values of z the first term in rhs(eq(z)) dominates. From this you get the idea to use the following instead:

fsolve(eq(z), H=h*sqrt(z^3*m))

@goli For large values of z the first term in rhs(eq(z)) dominates. From this you get the idea to use the following instead:

fsolve(eq(z), H=h*sqrt(z^3*m))

You have the variable z in series(gen, z, 8), but no variable z in gen, and I don't get the error message you get, but the result FAIL.

@Dennis_B If you solve numerically the easiest way to find M and integrals of M is to include the equation

eqM:= M(t)=z1(t)^2 * sin( z5(t) * t)^2;

#in the call to dsolve/numeric.

#like this

num:=dsolve({sys,eqM,z1(0)=z10,z2(0)=z20,z3(0)=z30,z4(0)=z40,z5(0)=z50},numeric,parameters=[A,d,g,k,m,v,z10,z20,z30,z40,z50],output=listprocedure);
#Choosing values for constants.
par:=[A = 1, d = 1, g = 10, k = 1, m = 1, v = 1, z10 = 0, z20 = 1, z30 = 0, z40 = 1, z50 = 0];
num(parameters=par);
Mn:=subs(num,M(t));
evalf(Int(Mn,0..5));
plot(Mn,0..10);

@Dennis_B If you solve numerically the easiest way to find M and integrals of M is to include the equation

eqM:= M(t)=z1(t)^2 * sin( z5(t) * t)^2;

#in the call to dsolve/numeric.

#like this

num:=dsolve({sys,eqM,z1(0)=z10,z2(0)=z20,z3(0)=z30,z4(0)=z40,z5(0)=z50},numeric,parameters=[A,d,g,k,m,v,z10,z20,z30,z40,z50],output=listprocedure);
#Choosing values for constants.
par:=[A = 1, d = 1, g = 10, k = 1, m = 1, v = 1, z10 = 0, z20 = 1, z30 = 0, z40 = 1, z50 = 0];
num(parameters=par);
Mn:=subs(num,M(t));
evalf(Int(Mn,0..5));
plot(Mn,0..10);

@smokeybob A particular solution is usually defined as any solution to the inhomogeneous equation. Thus replacing the arbitrary constants by any values (e.g. zero) will give you a particular solution.

@smokeybob A particular solution is usually defined as any solution to the inhomogeneous equation. Thus replacing the arbitrary constants by any values (e.g. zero) will give you a particular solution.

@Peter_Parker Sorry no typing error,  but the order in series is the problem. See my answer to the question asked separately.

@Peter_Parker Sorry no typing error,  but the order in series is the problem. See my answer to the question asked separately.

I was puzzled by the empty list in the procedure

proc(x,y) []; evalf(y+I*10^(-x)); end proc

What purpose does that serve?