It seems to me that a lot of the people reporting this problem are Danes (so am I). Is that a coincidence or what?

I haven't had that kind of problem for a very long time, maybe because I (almost) always remove all output before saving.

It seems to me that a lot of the people reporting this problem are Danes (so am I). Is that a coincidence or what?

I haven't had that kind of problem for a very long time, maybe because I (almost) always remove all output before saving.

LeastSquares work if B is converted to rational.

Bex:=convert(B,rational);
GaussianElimination(<A|Bex>);
Yp:=LeastSquares(A, Bex,free=t);
Norm(A.Yp-Bex,2);
Y:=LeastSquares(A, Bex,optimize);
Norm(A.Y-Bex,2);
evalf(%);

#Compare with your result from DirectSearch:

Yds:=<.899555775576182, .897427164193849, .896594849107632, .899558830317534>;
Norm(A.Yds-B,2);

LeastSquares work if B is converted to rational.

Bex:=convert(B,rational);
GaussianElimination(<A|Bex>);
Yp:=LeastSquares(A, Bex,free=t);
Norm(A.Yp-Bex,2);
Y:=LeastSquares(A, Bex,optimize);
Norm(A.Y-Bex,2);
evalf(%);

#Compare with your result from DirectSearch:

Yds:=<.899555775576182, .897427164193849, .896594849107632, .899558830317534>;
Norm(A.Yds-B,2);

@k4walker You are right. However, I did get 0 when I tried, but that must have been due to the fact that the order of the solutions is not guaranteed to be the same each time. So now I tried:

myDeterminant := (1/2)*l[4]*l[1]*(l[4]*cos(x[2]-2*x[3])+l[4]*cos(x[2])-l[1]*sin(2*x[2]-x[3])-l[1]*sin(x[3]));
solutions := solve(myDeterminant = 0);
eval~(myDeterminant,[solutions]);
simplify(%);
expand(%);
simplify(%);

From this we see that 5 of the 6 solutions found are correct.

@k4walker You are right. However, I did get 0 when I tried, but that must have been due to the fact that the order of the solutions is not guaranteed to be the same each time. So now I tried:

myDeterminant := (1/2)*l[4]*l[1]*(l[4]*cos(x[2]-2*x[3])+l[4]*cos(x[2])-l[1]*sin(2*x[2]-x[3])-l[1]*sin(x[3]));
solutions := solve(myDeterminant = 0);
eval~(myDeterminant,[solutions]);
simplify(%);
expand(%);
simplify(%);

From this we see that 5 of the 6 solutions found are correct.

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