@kh2n I misunderstood your question. Now if I understand you correctly, you just have to check if the answer provided by Maple in 'sol' is in fact correct, i.e. satisfies the ode and the bcs. This can be done by odetest:

odetest(sol,{ode, bcs});

The answer is {0}, which means that odetest says that ode and bcs are satisfied.

@kh2n I misunderstood your question. Now if I understand you correctly, you just have to check if the answer provided by Maple in 'sol' is in fact correct, i.e. satisfies the ode and the bcs. This can be done by odetest:

odetest(sol,{ode, bcs});

The answer is {0}, which means that odetest says that ode and bcs are satisfied.

@Red Horse No, the greatest in absolute value can be seen on the plot of the absolute value:

  plot(abs(diff(f,x,x)),x=-Pi/2..3*Pi/2);

@Red Horse No, the greatest in absolute value can be seen on the plot of the absolute value:

  plot(abs(diff(f,x,x)),x=-Pi/2..3*Pi/2);

You must have meant

E, V:= Eigenvectors(A);

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);