Did you mean

bcs1 := y(0) = 0, y(1) = 1, D(y)(1) = 0;

Notice two changes.

If so there is no problem.

I'm assuming that your ode is of first order since otherwise dsolve provides the answer directly.

There are a couple of ways of doing it. I begin with the one that comes to mind immediately, but which is not the simplest.

Isolate y'(x) from your equation. Then the other side is y'(x) in terms of y(x), which is known.
Here is an example, where y'(x) is already isolated.

eq:=diff(y(x),x)=sin(3*x+y(x))+cos(x^2)*y(x);
ic:=y(0)=1;
sol:=dsolve({eq,ic},numeric,output=listprocedure);
Y:=subs(sol,y(x));
Yp:=unapply(subs(y(x)=Y(x),rhs(eq)),x);
#Test:
Yp(1.2345);
#Two ways of plotting Y and Yp:
plot([Y,Yp],0..10);
plot([Y(x),Yp(x)],x=0..10);
#Without using Y or Yp and with or without 'output=listprocedure', odeplot can do the plotting:
plots:-odeplot(sol,[[x,y(x)],[x,rhs(eq)]],0..10,numpoints=400);
#Instead of using the numpoints option you can use the refine option (refine=1) if you include a range option (like range=0..10) in your call to dsolve/numeric.
####
#The simplest is to let dsolve do the construction of Yp:

sol2:=dsolve({eq,yp(x)=rhs(eq),ic},numeric,range=0..10);
plots:-odeplot(sol2,[[x,y(x)],[x,yp(x)]],0..10,refine=1);
#or without the need for isolating y'(x):

sol3:=dsolve({eq,yp(x)=diff(y(x),x),ic},numeric,range=0..10);

Maple has no problem solving the problem directly:

restart;
eq:=diff(F(x),x,x,x)+F(x)*diff(F(x),x,x)=0;
ibc:= F(-5)=0,D(F)(-5)=0,D(F)(5)=0.14;
res:=dsolve({eq,ibc},numeric,output=listprocedure);
plots:-odeplot(res,[x,F(x)],-5..5);
p1:=%:
##However, if you want alpha and beta in your splitting of the problem you can get them from the above result:

F0,F1:=op(subs(res,[F(x),diff(F(x),x)]));
F0(0);
F1(0);
#You can check that you get the correct solution by piecing g and h together:
gres:=dsolve({eq,F(-5)=0,D(F)(-5)=0,F(0)=F0(0)},numeric,output=listprocedure);
hres:=dsolve({eq,F(0)=F0(0),D(F)(0)=F1(0),D(F)(5)=0.14},numeric,output=listprocedure);
p2:=plots:-odeplot(gres,[x,F(x)],-5..0,color=blue):
p3:=plots:-odeplot(hres,[x,F(x)],0..5,color=blue):
plots:-display(p2,p3);
plots:-display(p1,p2,p3);

It seems more natural to me to collect all the data in a Matrix which you then export after the loops.
If that is acceptable, it can be done like this.

Initialize the result matrix Res and the row number r outside the loop:

Res:=Matrix():
r:=0:
#Then make the following changes:

 #dat:=Array([N,h1,occ]):
  dat:=Vector[row]([N,h1,occ]);
  #if occ<>0 then print(dat); fi:
   if occ<>0 then r:=r+1; Res(r,1..3):=dat fi:

Now you should be able to export the matrix Res.

A comment: For clarity you may as well put the procedure A1 before all the loops. However, none of its formal variables b1,e1,h0, h11, h12 is mentioned in the body! So the variables b, e, h, h1, h2 appearing in A1 refer to the global variables.
A better idea would be to avoid using a procedure A1 and do something like this:

restart:
with(LinearAlgebra): with(ExcelTools):
n:=0.8:    e:=1:   h2:=0.01: Res:=Matrix():
r:=0:
for N from 20 by 10 to 100 do
   for b from 0.01 by 0.05 to 0.26 do
      for h from 0.02 by 0.02 to 0.2 do
        B1:=Matrix(N,(i,j)->if i=j then e elif abs(i-j)=1 then h else 0 fi,datatype=float);
        B1[N,N-1]:=h2;
        B1[N-1,N]:=h2;
        B1[N,N]:=e-b;
        for h1 from 0.015 by 0.015 to 0.105 do
            B1[1,2]:=h1;
            B1[2,1]:=h1;
            p2:=Eigenvectors(B1,output='list'):
###And then complete as before ending with          
###
 .....
  dat:=Vector[row]([N,h1,occ]);
  if occ<>0 then r:=r+1; Res(r,1..3):=dat fi:
od: od: od: od:

You wouldn't have the problem if you used Maple input, i.e. 1D-input.

For 2D-input I have no I idea, as I don't use it.

Try replacing your ode by

ode2:=diff(x(t),t,t)=piecewise(x(t)=0,t,(t*sin(x(t))-x(t)^3)/sin(x(t)));

Notice that I use 'Array' instead of 'array'. (You have a misprint: Your first equation is eq1).

A:=dsolve({eq1,eq2,bcs},numeric, output=Array([0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,10]));
#The matrix containing all the numeric data is
M:=A[2,1];
#The matrix containing the numeric data for (eta, f(eta) ) is
M[1..,1..2];
#The matrix containing the numeric data for (eta, diff(f(eta),eta) ) is
M[1..,[1,3]];
#Here the numeric data for the two mentioned examples are exported:
ExportMatrix("F:/MapleDiverse/MaplePrimes/testF.txt",M[1..,1..2]);
ExportMatrix("F:/MapleDiverse/MaplePrimes/testFp.txt",M[1..,[1,3]]);

#The explanation if you need it:
T:=A[1,1];
T[1..2];
T[[1,3]];
#If you want to include the heading [eta f(eta) ] then instead export this matrix
<T[1..2],M[1..,1..2]>;

How about

newXfull:=Matrix(nodes,(i,j)->Xfull[i,j]/Xfull[1,j]);

or to do it your way:

newXfull:=Matrix([seq(zip(`/`,Column(Xfull,i),Row(Xfull,1)(i)),i=1..nodes)]);

Instead of making one animation you could make two. The first produced by a procedure F1 which doesn't plot the black points in the screens, the second produced by a procedure F2 which only plots those black points.
Then do

animate(F1, [t, t, t], t = 0 .. 60, frames = 100, scaling = constrained, background = Cameras, axes = frame); animate(F2, [t, t, t], t = 0 .. 60, frames = 100, scaling = constrained, axes = frame, trace = 100);
display(%%, %);

To get the data points you could save the results from the second animation:

res2 := animate(F2, [t, t, t], t = 0 .. 60, frames = 100, scaling = constrained, axes = frame, trace = 100):
A:=op([1,100],res2):
#The points in both screens:
L:=map2(op,1,A)[1..-2]:
nops(L);
#The points from the two screens separately:
L1:=[seq(L[2*i-1],i=1..100)]:
L2:=[seq(L[2*i],i=1..100)]:
display([seq(pointplot3d(L2[1..i]),i=1..100)],insequence=true,axes=boxed);
display([seq(pointplot3d(L[1..i]),i=1..200)],insequence=true,axes=boxed);

Is a real and between 0 and 1? Is N a positive integer? If so you can do:
restart;
Gamma := sqrt(1-a^2)/(1-a/z)*((1/z-a)/(1-a/z))^(N-1);
A:=Int(subs(z=exp(I*w),Gamma)*conjugate(subs(z=exp(I*w),Gamma)),w=-Pi..Pi);
As:=simplify(A) assuming a>0,a<1,N::posint;
value(As) assuming a>0,a<1;

restart;
ode := diff(theta(eta),eta,eta)+beta*diff(theta(eta),eta) = 0, theta(0) = 1, theta(5) = 0;
beta:=2;
N:=1:
#I used Array instead of array:
res:=dsolve({ode},numeric, output=Array([seq(i*0.1, i=0..10*N)]));
M:=res[2,1];
E,T,DT:=LinearAlgebra:-Column(M,1..3);
ET:=convert(zip(`[]`,E,T),list);
EDT:=convert(zip(`[]`,E,DT),list);
with(CurveFitting):
theta0:=PolynomialInterpolation(ET,eta);
dtheta0:=PolynomialInterpolation(EDT,eta);
plot(theta0,eta=0..2);
#But why use polynomial interpolation? Why not let dsolve/numeric do the whole job with output=listprocedure as follows:
resLP:=dsolve({ode},numeric, output=listprocedure);
#The procedures for theta and theta prime:
Th,DTh:=op(subs(resLP,[theta(eta),diff(theta(eta),eta)]));
Th(1.23);
DTh(1.23);

Both of these work:

readdata("F:/test.txt",[string,integer,integer,integer]);
evalindets(%,string,()->NULL);
#ImportMatrix may be easier to use:
ImportMatrix("F:/test.txt",delimiter=" ");
convert(%,listlist);
evalindets(%,name,()->NULL);

Since the "correct" answer is supposed to be 3 lists M1, M2, and M3, you really have 9 real variables.

It seems to me that you have 3 equations that are sums of dotproducts, the first one being

rho2.M1 - rho1.M2 = 0.

Now working with vectors and matrices instead of lists you could do as follows (if I understand your intention):

restart;
F := [x1/(1+t^2)+x2, -(x1+x2)*(1+1/(1+t^2)), -x3];
#Let the rho's be row vectors:
rho1 := <diff(F[1],x1)|diff(F[2],x1)|diff(F[3],x1)>;
rho2 := <diff(F[1],x2)|diff(F[2],x2)|diff(F[3],x2)>;
rho3 := <diff(F[1],x3)|diff(F[2],x3)|diff(F[3],x3)>;
theta := Matrix([[rho2, -rho1, <0|0|0>],[rho3, <0|0|0>, -rho1],[<0|0|0>, rho3,-rho2]]);
#The proposed solution is checked:
Mx:=<-1,-1,0,1,0,0,0,0,-1>; #The 3 solutions stacked on top of each other.
theta.Mx;
#Well, indeed a solution, but there are infinitely many more:
LinearAlgebra:-LinearSolve(theta,<0,0,0>,free=f);

I tried with success

restart;
evalf(Sum(Zeta(p-2.)*2.9^p/GAMMA(.5*p+1.), p = 1000..infinity));
evalf(Sum(Zeta(p-2.)*2.9^p/GAMMA(.5*p+1.), p = 4 .. 1000));
restart;
evalf(Sum(Zeta(p-2.)*8^p/GAMMA(.5*p+1.), p = 1000..infinity));
evalf(Sum(Zeta(p-2.)*8^p/GAMMA(.5*p+1.), p = 4 .. 1000));

But don't know why Maple doesn't succeed with the entire range at once.
#Comment: Actually the following works:
restart;
evalf[250](Sum(Zeta(p-2.)*2.9^p/GAMMA(.5*p+1.), p = 4 .. infinity));
evalf(%);
#But this doesn't. It returns Float(infinity) (and takes a long time)
evalf[2000](Sum(Zeta(p-2)*8^p/GAMMA(.5*p+1.), p = 4 .. infinity));

convert(par,listlist);
map( ((x, y)-> x^2+y^2)@op , %);