Clearly the answer from RealDomain:-solve({y=y,y<>0},y) is wrong and reveals a bug in RealDomain:-solve.

I think it can be traced from RealDomain:-solve via
SolveTools:-PolynomialSystemSolvers:-RealDomainFrontEnd
and
SolveTools:-PolynomialSystemSolvers:-RCRealSolver
to
SolveTools:-SemiAlgebraic,
which is being used here like this:
SolveTools:-SemiAlgebraic([0, y],{y});
It probably should have received the input [0,y<>0],{y} instead.
Indeed

SolveTools:-SemiAlgebraic([0, y<>0],{y});

returns [[y < 0], [0 < y]], which in the real domain is equivalent to y<>0.

@9009134 For your reformed program the system sys obtained by
 

sys,bcs:=selectremove(has,dsys3,diff);

still cannot immediately be solve for the highest derivatives hdiff:={diff(theta(x), x, x), diff(u(x), x, x, x, x), diff(w(x), x, x, x, x)}.
But sys[1] can be written
diff(u(x), x, x)+diff(w(x), x)*diff(w(x), x, x) = 0.
If that equation is replaced by the same differentiated twice the resulting new system can be solved for the highest deivatives hdiff.
Since you differentiated sys[1] twice you have 2 constraints:

diff(u(x), x, x)+diff(w(x), x)*diff(w(x), x, x) = 0 and
diff(u(x), x, x,x)+diff(w(x), x,x)^2 + diff(w(x), x)*diff(w(x), x, x,x) = 0.

Since the lhs of these will be constant you just need to evaluate them at one of the endpoints, e.g. 0
This gives you two new boundary conditions besides the 10 you already got. Thus you must remove two of the old ones.

Actually solving the new system for hdiff reveals that the denominator
(-1+400000*w(x))^6*(2*diff(u(x), x)+diff(w(x), x)^2)
appears.
This means that D(u)(0) and D(w)(0) cannot both be zero if you want to avoid a singularity at x=0. It's the same at x=0.001.
Thus in revising your bcs you should have this in mind.


 

A solution is indeed found in Maple 2016.2. It takes a while and the result is so enormous that it is not printed, but a message is printed:
`[Length of output exceeds limit of 1000000]`

The result is there, but you should be happy that it is not printed.
Suppose you save the solution in sol.
Then try e.g.
 

sol:=solve({e1, e2, e3, e4, e5}, {S__1, S__2, S__3, S__4, S__5});
nops([sol]); # 2
sol[1]; #Very simple
length(sol[2]); # 6216065
## Setting the parameters to 1 all of them for ease of testing:
params:=indets({e1, e2, e3, e4, e5},name) minus {S__1, S__2, S__3, S__4, S__5};
evalf[50](eval(sol[2],params=~1));

Actually it is the second version that is correct. The first is false.
 

restart;
eqns:=
[x=1, -1/(exp(h)+1/exp(h))^(1/2)/(exp(h)-1/exp(h))^(1/2)*(exp(h)^2-2+1/exp(h)^2)^(1/2)*x = tanh(h)^(1/2)];
eval(eqns,{x=1,h=2});
evalf(%); # Not equal

For h>0 eqns has no solutions for x. For h<=0 it does. (And nobody says that h couldn't be complex either).
Looking a little closer:
 

eq2:=eval(eqns[2],x=1);
rhs(eq2)^2-lhs(eq2)^2;
convert(%,exp); # 0

It may be that solve in getting the first incorrect result squares eq2 as I just did.

To get Vector/Matrix/Array output from dsolve/numeric use output=Array(...).
Example:
 

restart;
sys:=diff(x(t),t)=y(t),diff(y(t),t)+x(t)=sin(t);
ics:=x(0)=0,y(0)=1;
sol:=dsolve({sys,ics});
b:=10: N:=64: h:=b/(N-1);
T:=Array(1..N,i->i*h,datatype=float);
res:=dsolve({sys,ics},numeric,output=T);
A:=res[2,1];
plot(A[..,1..2]); # Test
X:=A[..,2]; # X values as a vector
Y:=A[..,3]; # Y values as a vector
plot(T,Y); # Another test

For the boundary conditions you have Maple's dsolve/numeric/bvp has severe problems and that may be because the problem has no solution. For the conventional boundary values f(0)=0, D(f)(0)=0, D(f)(inf) = k (k some constant, e.g. 2) there is no problem:
 

restart;
ode:=diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)=0
bcs1:=D(f)(0)=1,(D@@2)(f)(0)=0,(D@@2)(f)(6)=2; #Your boundary conditions
dsolve({ode,bcs1},numeric,initmesh=8192); # No success even with initmesh high
## Conventional boundary conditions:
bcs2:=f(0)=0,D(f)(0)=0,D(f)(6)=2;
res2:=dsolve({ode,bcs2},numeric);
plots:-odeplot(res2,[seq([eta,diff(f(eta),[eta$i])],i=0..2)]);

Since dsolve/numeric/bvp has problems I'm not at all surprised that your scheme does too. On the contrary that should be expected.
To examine your scheme you may try using the boundary conditions bcs2 instead. If your scheme doesn't work quite as then look for logical mistakes.

The difference between the 'good' and the 'bad' solves the corresponding homogeneous pde.
Both solutions are correct.
 

pdetest(good,eqn_sys); # OK
pdetest(bad,eqn_sys); # OK
##
solhom:=y(t,x)=op(combine(expand(rhs~(good-~bad))));
HOM:=D[1, 1](y)(t, x)+y(t, x) = 0;
pdetest(solhom,HOM); # OK

Although Kitonum has already answered your question perfectly fine I shall add mine too.

Your question may be related to your question of June 19 2018:
https://mapleprimes.com/questions/224941-Error-in-Dsolvenumericbvpconvertsys#answer250475 

Now you give us 3 odes with only 2 unknown functions. Two of the odes are linear with constant coefficients thus basically trivial:
 

ode1:=1.857142857*10^11*diff(u(x), x, x)=0;
ode2:=1.547619048*10^10*diff(w(x), x, x, x)-1.300000000*10^11*diff(w(x), x)=0;

In addition we have the homogeneous boundary conditions:
 

bcs:= u(0)=0, u(0.001)=0, w(0)=0, w(0.001)=0, D(w)(0)=0;

Using dsolve on that system:

dsolve({ode1,ode2,bcs});

gives you (not surprisingly) u(x)=0, w(x)=0.
But you have a third ode (call it ode3) that only involves w. That one doesn't have the trivial solution w(x) = 0 for all x:

eval(ode3,w=0);

gives you .3693149535 = 0.


 

Copying the examples from the CUDA help page I get at the given stage in Maple 2018.1 as well as in Maple 2017.3:

tCUDA := time[real](M1.M2);
Error, (in unknown) CUBLAS internal error

In Maple 2016.2, 2015.2, 18,  17, 16, and 15 there is no problem with the example section for CUDA on my machine. Those releases of Maple are all installed on the same computer.

I shall submit an SCR.

Although as I pointed out in a comment, x>a and x<>c evaluates to x>a the equivalent formulation
 

x>a and (x>c or x<c);

does not. It returns unevaluated.
And so does:
 

x>-infinity and x<infinity and (x>1/2 or x<1/2);

solve handles the latter just fine:

solve(x>-infinity and x<infinity and (x>1/2 or x<1/2),x);

and returns RealRange(Open(1/2), Open(infinity)), RealRange(Open(-infinity), Open(1/2)).

In general x<>c is not equivalent to (x<c or x>c) since the latter assumes an order where the first simply means 'x is different from c'.

As mehdi jafari is saying you should consider the "assignment" in the loop an equation. Call that equation eq[i,j], say.
Forget about doing a loop.
Make the set of equations S:={seq(seq( eq[i,j], i=1..N),j=0..N)};
Use solve to find the unknown u[i,j].
Give us the code as text or upload a worksheet if you need further assistance.

Since you could solve the problem exactly, dsolve ought to be able to do it too.
Following your lead dsolve can do it stepwise using that the first ode only has f[1], the second only f[1] and f[2], etc.
 

restart;
N:=4; alpha:=5*3.14/180; r:=10; Ha:=5; H:=1;
Rf:=diff(f[m-1](x),x,x,x)+2*alpha*r*sum(f[m-1-n](x)*diff(f[n](x),x),n=0..m-1)
+(4-Ha)*(alpha)^2*diff(f[m-1](x),x); 
ode:=diff(f[m](x),x,x,x)-CHI*diff(f[m-1](x),x,x,x)=h*H*Rf;
CHI:=`if`(m>1,1,0);
f[0]:=x->1-x^2;
## Using a list to keep the order:
sysL:=[seq(ode,m=1..N)];
bcsL:=[seq({f[m](0)=0,D(f[m])(0)=0,f[m](1)=0},m=1..N)];
for m from 1 to N do
  res[m]:=dsolve(eval({sysL[m]},{seq(res[i],i=1..m-1)}) union bcsL[m])
end do;
F:=f[0](x)+add(rhs(res[m]),m=1..N);
plot(eval(F,h=-0.9),x=0..1);
plot(eval(diff(F,x),x=1),h=-1.5..-0.2);

The plots look like the plots done with the numerical solution.

I tried your system and ended up using numerical solution.
 

restart;
N:=4; alpha:=5*3.14/180; r:=10; Ha:=5; H:=1;
Rf:=diff(f[m-1](x),x,x,x)+2*alpha*r*sum(f[m-1-n](x)*diff(f[n](x),x),n=0..m-1)
+(4-Ha)*(alpha)^2*diff(f[m-1](x),x); 
####
ode:=diff(f[m](x),x,x,x)-CHI*diff(f[m-1](x),x,x,x)=h*H*Rf;
CHI:=`if`(m>1,1,0);
f[0]:=x->1-x^2;
#### The ode system for f[i], i=1..N
sys:={seq(ode,m=1..N)};
bcs:=`union`(seq({f[m](0)=0,D(f[m])(0)=0,f[m](1)=0},m=1..N));
#### Introducing F as the sum avoiding using f (avoiding the possible confusion with f[i])
eq:=F(x)=add(f[m](x),m=0..N);
#### Since D(F)(1) is needed we introduce a second order ode for F:
odeF:=diff(eq,x,x);
#### Initial conditions for F:
icF0:=eval[recurse](eq,{x=0} union bcs);
eq1:=convert(diff(eq,x),D);
icF1:=eval[recurse](eq1,{x=0} union bcs);
#### The total system
SYS:=sys union {odeF}:
BCS:=bcs union {icF0,icF1}:
#### To obtain the two plots we let res return a list of two results:
res:=proc(h1) local sol,F1;
  if not h1::realcons then return 'procname(_passed)' end if;
  sol:=dsolve(eval(SYS,h=h1) union BCS,numeric,output=listprocedure);
  F1:=subs(sol,diff(F(x),x))(1);
  [sol,F1]
end proc;
#### Testing
res(-0.9);
#### The plots:
plots:-odeplot(res(-0.9)[1],[x,F(x)],0..1);
plot(res(h)[2],h=-1.5..-0.2,caption="D(F)(1) as a function of h");

The attempt to solve sys with bcs exactly fails:
 

dsolve(sys union bcs); #No answer

This happens if you use 2D input but doesn't if you use 1D input (a.k.a. Maple input).
I tried both in Maple 2018.1 and in Maple 2017.3.
No problem when using 1D input, but with 2D input I got the error message

Error, invalid loop statement termination

Yet another irritating difference between 1D input and 2D input.

MaplePrimes18-07-11_2D_problem.mw

Note. Just to make sure the 2D input problem has nothing to do with the fact that the statement doesn't do anything I tried this version:
 

if true then
  for i from 1 to 2 do
    i
  end do
else
  6
end if;

The same conclusion: Fine in 1D, error in 2D.

For numeric solution your system cannot contain vectors. You must rewrite the vector system so you have a system of scalar equations. That ought to be easy.
For help on this you must supply either code in text form or a worksheet (upload using the fat green arrow in the MaplePrimes editor).