This is clearly difficult to deal with.
I changed blt from 15.5 to just 3. I had only success with s = 0.
It appears that theta for all practical purposes is zero for eta > 0.003. Thus I tried solving Eq1 with theta = 0.
That was very easy. After that my intention was to use the results for f to solve Eq2.

restart;
Digits:=15:
with(plots):
#blt:=15.5:
A:=1:
M:=1:
lambda:=1:
fw:=1:
rhos:=8933:
rhof:=997.1:
rhob:=14918.11:
rhob2:=20939.1:
kn:=0.613:
kf:=401:
cps:=385:
cpf:=4179:
Pr:=6.2:
Eq1 := (1/(1-s)^(2.5))*diff(f(eta),eta,eta,eta)-(1-s+s*(rhos/rhof))*(diff(f(eta),eta)^(2)-f(eta)*diff(f(eta),eta,eta)+A*(diff(f(eta),eta)+(1/2)*eta*diff(f(eta),eta,eta)))-M*diff(f(eta),eta)+(1-s+s*((rhob/rhob2)))*lambda*theta(eta)=0;
Eq2 :=(1/Pr)*((kn/kf)/(1-s+s*(rhos*cps/rhof*cpf)))*diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-diff(f(eta),eta)*theta(eta)-A*(2*theta(eta)+0.5*eta*diff(theta(eta),eta))=0;
blt:=3:
bcs1 := f(0) = fw, (D(f))(0) = 1, (D(f))(blt) = 0, theta(0) = 1, theta(blt) = 0;
L := [0,0.05,0.1,0.2];
R:=dsolve(eval({Eq1, Eq2, bcs1}, s = 0),numeric, output = listprocedure,initmesh=1024,maxmesh=8192);
odeplot(R,[eta,theta(eta)],0..0.003);
odeplot(R,[[eta,f(eta)],[eta,theta(eta)]],0..blt);
#################
##Now the approximation obtained by setting theta = 0 in Eq1.
## I'm saving the procedures for f, f', f'' for use later:
blt:=8: #I could get away with blt=8
AP:=NULL:
for k from 1 to 4 do
  Rf[k] := dsolve(eval({Eq1f, f(0)=1,D(f)(0)=1,D(f)(blt)=0}, s = L[k]),numeric, output = listprocedure,AP);
  AP:=approxsoln=Rf[k];
  F[k],F1[k],F2[k]:=op(subs(Rf[k],[seq(diff(f(eta),[eta$i]),i=0..2)]))
  end do;
display(seq(odeplot(Rf[k],[eta,f(eta)],0..blt),k=1..4),insequence);
################
#Now I tried the following loop, but had success only with k=1, i.e. s=0:
AP:=NULL:
for k from 1 to 4 do
 Rt[k]:=dsolve({subs(diff(f(eta),eta)=F1[k](eta),f(eta)=F[k](eta),eval(Eq2,s=L[k])),theta(0)=1,theta(0.003)=0},numeric,known=[F[k],F1[k]],AP);
AP:=approxsoln=Rt[k]
end do;
odeplot(Rt[1],[eta,theta(eta)],0..0.003);
###############FINALLY I tried an approximation on Eq2 which uses that theta decreases steeply to zero.
## I'm using that Eq2 can be approximated by K*theta''+f*theta' = 0, but also by
## K*theta''+(f*theta)' = 0 since f '*theta is relatively small.
## Integrating once using that theta and theta' are zero for eta = 0.003 we get the first order ode:
Eq20 :=(1/Pr)*((kn/kf)/(1-s+s*(rhos*cps/rhof*cpf)))*diff(theta(eta),eta)+f(eta)*theta(eta)=0;
#This can be solve symbolically in terms of f:
solTH:=dsolve({Eq20,theta(0)=1},theta(eta));
#Now we can plot theta
for k from 1 to 4 do
  if k=1 then r:=0.003 else r:=1e-8 end if;
  pp[k]:=plot(eval(rhs(solTH),{s=L[k],f=F[k]}),eta=0..r,thickness=5)
end do:
display(Array([seq(pp[i],i=1..4)]));
###In fact we may as well replace f with 1 in solTH. It makes no difference of any importance to the results, but has the grat advantage that solTH then can be integrated without knowing f:
solTH2:=value(eval(solTH,f=1));
#Now on the other hand we could replace theta in Eq1 with that explicit expression:
Eq1f2:=subs(solTH2,Eq1);
##Now use a dsolve loop as above taking blt=5:
blt:=5:
AP:=NULL:
for k from 1 to 4 do
  Rf2[k] := dsolve(eval({Eq1f2, f(0)=1,D(f)(0)=1,D(f)(blt)=0}, s = L[k]),numeric, output = listprocedure,AP,maxmesh=8192,abserr=1e-5);
  AP:=approxsoln=Rf2[k];
  #F[k],F1[k],F2[k]:=op(subs(Rf2[k],[seq(diff(f(eta),[eta$i]),i=0..2)]))
  end do;
The graphs are really indistinguishable from the ones obtained using Eq1f with blt=8 (ignoring eta>5).

My conclusion would be to use solTH2 as the solutions for theta and as the solutions for f the ones obtained from Eq1f.

I don't think you are going to find a solution with the values of R and delta that you have. So either R has to be smaller or delta has to be larger.
Instead of going through a loop solving the ode system 3000 times I used the parameters option in dsolve and also fsolve.
For convenience I bring the entire code.

restart;
R:=0.3:
theta_min:=Pi/6:
theta_max:=Pi/2:
betta_max:=evalf(Pi/180*80);
p:=2*10^5:
theta0:=s->Pi/3/s_end*s+Pi/6:
r0:=s->R*sin(theta0(s)):
s_end:=evalf(R*(theta_max-theta_min)):
sol1:=solve({sin(betta_max)=c/r0(0)},{c});
const1:=0.1477211630;
betta0:=s->arcsin(const1/r0(s)):
betta:=s->arcsin(r(s)/r0(s)*sin(betta0(s))):
A:=s->cos(betta(s))/cos(betta0(s)):
T1:=s->rT1(s)/r(s):
T2:=s->T1(s)*tan(betta(s))^2:
delta:=0.001:
sys := diff(rT1(s), s)-A(s)*T2(s)*cos(theta(s)),diff(theta(s), s)-A(s)/T1(s)*(p-T2(s)*sin(theta(s)/r(s))),diff(r(s),s)-A(s)*cos(theta(s)),diff(z(s),s)-A(s)*sin(theta(s));
rT1_n:=p*Pi*r_min^2/2/Pi/sin(theta_min);
##### No changes made above to the lines included. ############
## Now use the parameters option with r_min as parameter:
F := dsolve({sys,rT1(0)=rT1_n, theta(0)=theta_min,r(0) = r_min, z(0) = 0}, numeric,output=listprocedure,parameters=[r_min]);
r_ravn:=subs(F,r(s)); #The procedure for evaluating r(s) extracted
#Testing. Setting the parameter:
F(parameters=[0.152]);
F(s_end); #Values at s_end
##Plotting (scaling so all 3 appear) (still r_min=0.152)
plots:-odeplot(F,[[s,10*r(s)],[s,rT1(s)/20000],[s,theta(s)]],0..s_end);
## Another plot with separate plots (still r_min=0.152)
use plots in display(Array([seq(odeplot(F,[s,fu(s)],0..s_end),fu=[r,rT1,theta])])) end use;
## Making a procedure which sets the parameter and plots. Default all 3 functions as an array.
q:=proc(r_min,L::list:=[r,rT1,theta]) uses plots; F(parameters=[r_min]);
  display(Array([seq(odeplot(F,[s,fu(s)],0..s_end),fu=L)])) end proc;
q(.001); #Plotting all with r_min = 0.001
q(0.1,[r]); #Plotting only r(s) with r_min=0.1
plots:-animate(q,[rmin],rmin=0.01..0.152,frames=100); #Animation (all 3)
plots:-animate(q,[rmin,[r]],rmin=0.001..0.152,frames=100); #Animation (only r)
##It appears that r stays below R-delta for all s = 0..s_end
## Thus I replace abs(r_ravn(s_end)-R)=delta by r_ravn(s_end)-R+delta=0 below.
## In the procedure pp the default is to use the existing delta, but another can be given as a second argument:
pp:=proc(r_min,dd::numeric:=delta) F(parameters=[r_min]); r_ravn(s_end)-R+dd end proc;
#Testing:
pp(.152);
pp(.152,.01);
pp(.152,.02);
#Now use fsolve. I expected no success with the default delta and didn't have any:
fsolve(pp,0.01..0.152); returns unevaluated
fsolve(rm->pp(rm,.02),0.01..0.152); #Success
### Finally to see that R-delta is too large we can plot r_ravn(s_end) as a funcion of r_min.
## The plot clearly shows that the maximal value for r_ravn(s_end) is for r_min=0.152:
vv:=proc(r_min) F(parameters=[r_min]); r_ravn(s_end) end proc;
plot(vv,0.01..0.152);
vv(0.152);

The second time around z is no more a symbol.
I understand this to be a small version of something more complicated. Is there a need to introduce symbols in Z?
###
If not then of course there is no problem:
restart;
for k from 1 to 2 do Z:=Matrix(2,2); Z[1,1]:=1; print(Z); end do;
###
Anyway to solve the problem as presented you could do either one of these:

restart;
for k from 1 to 2 do Z:=Matrix(2,2,symbol=z); z[1,1]:=1; print(Z); z:='z'; end do;
###
restart;
for k from 1 to 2 do Z:=Matrix(2,2,symbol=z); Z[1,1]:=1; print(Z); end do;

Do you actually have a nonlinear eigenvalue problem, where unfortunately for dsolve you have already found an eigenvalue and that is inserted into your system after which there is no parameter only concrete values for the constants.
If so it is better to let Maple find an eigenvalue and a corresponding solution at the same time.

Let me start by giving a very simple example after which I shall make an experiment with your problem.
restart;
ode:=diff(x(t),t,t)+lambda*x(t)^2=0; #Nonlinear ode
res:=dsolve({ode,x(0)=0,x(1)=0,D(x)(0)=1},numeric); #Extra requirement D(x)(0)=1 <> 0
res(0); #You see the eigenvalue
plots:-odeplot(res,[t,x(t)],0..1); #Corresponding solution for x
#Now insert the eigenvalue found into ode:
ode1:=subs(res(0)[-1],ode);
#Try solving with homogeneous boundary conditions:
res1:=dsolve({ode1,x(0)=0,x(1)=0},numeric);
res1(0.5); #Indicates clearly enough already that you just have the zero solution
res2:=dsolve({ode1,x(0)=0,D(x)(0)=1},numeric); #Using the extra condition instead of x(1)=0.
plots:-odeplot(res2,[t,x(t)],0..1);
################################
Now let us try an experiment with your concrete system containing no parameter lambda.
We shall rather arbitrarily insert a parameter lambda into your system and shall add an extra (nonzero) boundary condition.

I use the same ideas I used in my other answer, but I repeat it here.

restart;
dsys3 := {0.1873004296e-2*(diff(f3(x), x, x, x, x))-.3714109950*(diff(f3(x), x, x))+13.96314014*(diff(f1(x), x))+1.822500000*10^(-10)*(diff(f2(x), x, x))-4.888800000*10^(-7)*f2(x)+22.72337812*f3(x) = 0, -901.7765285*(diff(f1(x), x, x))+3.153274902*10^(-10)*(diff(f2(x), x, x, x))+0.4807756008e-5*(diff(f2(x), x))-27.92641797*(diff(f3(x), x))+18074.25057*f1(x)-3.100642108*(diff(f3(x), x, x))*(diff(f3(x), x))+7.74301535*(diff(f3(x), x))*f3(x) = 0, 0.2249991345e-1*(diff(f2(x), x, x, x, x))-227.1861332*(diff(f2(x), x, x))+29537.78948*f2(x)-2.841067252*10^(-11)*(diff(f1(x), x, x, x))-3.392672540*10^(-7)*(diff(f1(x), x))+7.650000000*10^(-11)*(diff(f3(x), x, x))-1.400000000*10^(-7)*f3(x)-5.500000000*10^(-7)*f3(x)^2+9.278333333*10^(-13)*(diff(f3(x), x))^2+3.300000000*10^(-13)*(diff(f3(x), x, x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, (D(f2))(0) = 0, (D(f2))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0};
for i to 3 do dsys3[i] end do;
diff(dsys3[2],x);
eq:=subs(solve(dsys3[3],{diff(f2(x),x$4)}),%);
sys:=solve({dsys3[1],dsys3[3],eq},{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
dsys3[4..-1];
nops(%);
cdn:=eval(convert(dsys3[2],D),x=0);
#The insertion of lambda will be done here (obviously this is completely arbitrary and just an illustration):
sysE:={sys[1],sys[2],lhs(sys[3])=lambda*rhs(sys[3])};
resE:=dsolve(sysE union dsys3[4..-1] union {cdn, D(f1)(0)=1},numeric);
resE(0);
#Eigenvalue is -.493863843115919 i.e. not 1 so we don't have your original system.
#However, I inserted this lambda arbitrarily without knowing your original problem, so lambda=1 was hardly to be expected.
plots:-odeplot(resE,[[x,f1(x)],[x,f2(x)],[x,f3(x)]],0..1,thickness=3);

Now that we have the constants, we can begin by observing that solving for the highest derivatives isn't possible.

solve(dsys3[1..3],{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
##However by manipulating the system some, e.g. as follows you can solve the resulting system for the highest derivatives.
Please pay attention to if the ordering in the set dsys3[1..3] is such that dsys3[2] is the ode without any 4th derivative. It is assumed in the following:
Since we are differentiating dsys3[2] there will be the constraint dsys3[2] to keep in mind.
restart;
dsys3 := {0.1873004296e-2*(diff(f3(x), x, x, x, x))-.3714109950*(diff(f3(x), x, x))+13.96314014*(diff(f1(x), x))+1.822500000*10^(-10)*(diff(f2(x), x, x))-4.888800000*10^(-7)*f2(x)+22.72337812*f3(x) = 0, -901.7765285*(diff(f1(x), x, x))+3.153274902*10^(-10)*(diff(f2(x), x, x, x))+0.4807756008e-5*(diff(f2(x), x))-27.92641797*(diff(f3(x), x))+18074.25057*f1(x)-3.100642108*(diff(f3(x), x, x))*(diff(f3(x), x))+7.74301535*(diff(f3(x), x))*f3(x) = 0, 0.2249991345e-1*(diff(f2(x), x, x, x, x))-227.1861332*(diff(f2(x), x, x))+29537.78948*f2(x)-2.841067252*10^(-11)*(diff(f1(x), x, x, x))-3.392672540*10^(-7)*(diff(f1(x), x))+7.650000000*10^(-11)*(diff(f3(x), x, x))-1.400000000*10^(-7)*f3(x)-5.500000000*10^(-7)*f3(x)^2+9.278333333*10^(-13)*(diff(f3(x), x))^2+3.300000000*10^(-13)*(diff(f3(x), x, x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, (D(f2))(0) = 0, (D(f2))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0};
for i to 3 do dsys3[i] end do; #JUST FOR HAVING A LOOK at the equations
diff(dsys3[2],x);
eq:=subs(solve(dsys3[3],{diff(f2(x),x$4)}),%);
sys:=solve({dsys3[1],dsys3[3],eq},{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
dsys3[4..-1]; #The boundary conditions
nops(%);
cdn:=eval(convert(dsys3[2],D),x=0); #The above mentioned constraint
res:=dsolve(sys union dsys3[4..-1] union {cdn},numeric);
res(0.5);
plots:-odeplot(res,[x,f3(x)],0..1);
#Now belatedly we notice that there is a trivial solution:
odetest({f1(x)=0,f2(x)=0,f3(x)=0},dsys3[1..3]);
############
Thus you are I assume looking for the possibility of a nontrivial solution?

I think the small parameter method is something like the following, where because of only two new independent variables x1 and x2 we go to second order in epsilon. It is important of course to notice that x(t) = 1 is a solution (in particular when epsilon=0).

restart;
eq:= diff(x(t),t$2) - epsilon*(alpha*x(t)^4 - beta)*diff(x(t),t) - x(t) + x(t)^3;
eq2:= eval(eq, x(t)= 1 + epsilon*x[1](t) + epsilon^2*x[2](t));
collect(eq2,epsilon,factor);
ode1:=coeff(%,epsilon);
ode2:=coeff(%%,epsilon^2);
res:=dsolve({ode1,ode2},explicit);
#To shorten:
combine(res);

#################Numerical comparison for epsilon = 0.1 ##########
res1:=dsolve({eval(eq,{alpha=1,beta=1,epsilon=.1}),x(0)=1.11,D(x)(0)=0},numeric); #1.11 = 1+eps+eps^2
plots:-odeplot(res1,[t,x(t)],0..15); p1:=%:
##Using numerical solution here for convenience in the comparison
res2:=dsolve({ode1,eval(ode2,{alpha=1,beta=1}),x[1](0)=1,x[2](0)=1,D(x[1])(0)=0,D(x[2])(0)=0},numeric);
plots:-odeplot(res2,[t,1+0.1*x[1](t)+0.01*x[2](t)],0..15,color=blue); p2:=%:
plots:-display(p1,p2);
##To use the symbolic solution for x1,x2 just remove 'numeric':
res2a:=dsolve({ode1,eval(ode2,{alpha=1,beta=1}),x[1](0)=1,x[2](0)=1,D(x[1])(0)=0,D(x[2])(0)=0});
X:=subs(res2a,1+0.1*x[1](t)+0.01*x[2](t));
p3:=plot(X,t=0..15,color=blue);
plots:-display(p1,p3);

Was your intention with p something like this:

simplify(piecewise(abs(z)>6^(2/3)/4,piecewise(z<0, eq30, eq31),piecewise(z<0, eqc20, eqc21)));
p:=unapply(%,z,t);

eq:=diff(y(x),x)=cos(x)+y(y(x)-2);
res:=dsolve({eq,y(0)=1},numeric,delaymax=1);
plots:-odeplot(res,[x,y(x)],0..1);

#If you try res(2) you get an error message saying:
res(2);
Error, (in res) cannot evaluate the solution further right of 1.0058110, probably a singularity

You could use Spline from the CurveFitting package:

res:=CurveFitting:-Spline(Data1,x);
m,M:=(min,max)(map2(op,1,Data1));
plot(Data1,style=point,symbolsize=20,color=blue); p1:=%:
plot(res,x=m..M); p2:=%:
plots:-display(p1,p2);

Somehow `dsolve/numeric/bvp/convertsys` fails to isolate YP[4] when converting to a first order system. I shall show that below. YP[4] stands for diff(f(eta),eta$4).
A workaround is to do this isolation before we try dsolve.

Thus e.g. one could do:

ode12:=solve({ode1,ode2},{diff(f(eta),eta$4),diff(theta(eta),eta,eta)});
sys:=eval(ode12 union {bcs},E);
res:=dsolve(eval(sys,epsilon=L[1]),numeric,method=bvp[midrich],maxmesh=1024);
plots:-odeplot(res,[eta,f(eta)],0..5);

##### Another workaround (simpler):
sys:=convert(sys,rational);
## A third workaround:
sys:=collect(sys,diff);

####################
To exhibit the problem you can begin with the redefined 'sys' from above:
`dsolve/numeric/bvp/convertsys`(eval(sys[1..2],epsilon=L[1]),sys[3..-1],[f(eta),theta(eta)],eta,cont);
#Obviously here YP[4] is isolated on the left, since it started out that way.
# Now try the nonisolated system:
sys:=eval([ode1, ode2, bcs], E);
`dsolve/numeric/bvp/convertsys`(eval(sys[1..2],epsilon=L[1]),sys[3..-1],[f(eta),theta(eta)],eta,cont);
## You will notice that for i=1,2,3,5,and 6, YP[i] has been isolated, but an equation defining YP[4] implicitly has been left unsolved. I don't believe that that is intentional: I tried changing f(0)=0 to f(0)=1 and that didn't make any difference to convertsys.
So my guess is that it is indeed a bug in `dsolve/numeric/bvp/convertsys`.
####
It appears that the problem is in line 55 of `dsolve/numeric/bvp/convertsys` where a (global) procedure `dsn/solve` fails to isolate YP[4].
Finally and ironically, the problem in `dsn/solve` is with the isolate command in line 30 where isolate cannot isolate YP[4].

Under the assumption that b^2-4*c < 0 clearly the integral must be real (except of course for an additive constant: the integral is indefinite).
As an example I plot the value for b=2, c=2, where b^2-4*c = -4 < 0:
res:=int((y^4+b*y^2+c)^(-1/2),y);
plot(eval(res,{b=2,c=2}),y=-8..8);

Notice that the second derivatives of alpha and beta only appear in eq2. Thus solving (algebraically) for the second derivatives cannot be done as in

solve({eq1,eq2}, diff({alpha(t),beta(t)},t,t));

By using the word 'algebraically' I mean that if we consider the derivatives of alpha and beta of orders 0,1,2 as just six variables (like a0,a1,a2,b0,b1,b2) we cannot solve for {a2,b2} in terms of {a0,a1,b0,b1}.
We can get around it though: See the second method below.

Thus in order to find expressions for the second derivatives we could differentiate eq1 or let Maple do it. However, in differentiating eq1 information is lost as any constant on the right becomes a zero. Thus this equation must be kept at least evaluated at one point.

Anyway, Maple can do this itself like this:

solution := dsolve([eq1, eq2, CI], numeric, range = 0 .. 2,method=rkf45_dae);
#You will get an error about the initial conditions not satisfying the constraint, and in fact they don't:

eval(convert(eq1,D),t=0);
eval(%,{CI});

Thus you must clear up that problem, i.e. change the initial conditions.
################################ Alternative solution ################
Alternatively you can do like the following. Here we avoid any differentiation.
Notice that no requirements can be put on the derivative of beta at zero. This reflects the constraint problem above.
eq_beta:=solve(eq1,{diff(beta(t),t)});
eq3:=eval(eq2,%);
solution:=dsolve({eq3,op(eq_beta),alpha(0) = 0, beta(0) = 0,D(alpha)(0) = w[0]},numeric);
plots:-odeplot(solution,[[t,alpha(t)],[t,beta(t)]],0..0.9);

The _Z in the RootOf result is just a dummy variable, the unknown in the expression inside RootOf. The meaning is that the result is a root in that expression, but Maple either could not express it analytically or didn't 'want' to unless forced to as in the case of the 4 solutions to a quartic polynomial (see allvalues).
In your case there is no chance for an analytical solution (i.e. a  formula for the solution).

If you set

res:=solve(equ,c);

then you can use evalf to get a numerical evaluation of (one of) the root(s):
evalf(res);

You will find a result that looks like it is likely that c=3100 is an exact solution. It is clear from your equation that c=3100 is indeed a solution, since 3100^2 = 9610000, so c=31000 makes both sides of the equation zero.

This solution is apparently not the one you want, so you could try fsolve with a start guess for c:

fsolve(equ,c=13000);
You will find a c-value close to zero. Not the one you want either.

For the square roots all to be real you need  c > 6300 thus it makes sense to simplify equ to:

equa:=simplify(equ) assuming c>6300;

You could try plotting the difference between the two sides (of equ or equa):
plot((rhs-lhs)(equa),c=6300..15000);

After that you should ask yourself if there are indeed any solutions or if you just wrote down the equation incorrectly.

The operator Mf of multiplication by the function f can be defined like this:

Mf:=g->f*g;
#Applying Mf to a function h (in the mathematical sense, not Maple)
(Mf(h))(x); #The parenthesis around Mf(h) i.e. (Mf(h)) are superfluous, but are added for clarity of meaning.
((D-Mf)@@2)(g)(x);
expand(%);

#An operator M which maps functions onto multiplication operators:

M:=f->(g->f*g);
#Thus M(f) is the operator which acts like this:
M(h)(g)(x);
((D-M(f))@@2)(g)(x);
expand(%);

########Multiplying by a constant function:
M(89)(g)(x);
a:=x->k; #The function a defined by a(x)=k for all x
M(a)(g)(x);

############## Followup question
#To accomplish your intended result you can use M above:

restart;
M:=f->(g->f*g);
id:=(x,m)->x;
M(id)(g)(y,n); #test
(M(id)@D[1])(g)(x,m); #test
((M(id)@D[1])@@2)(F)(x,m);

#Alternatively use id but not M:
Mx:=g->id*g;
((Mx@D[1])@@2)(F)(x,m);

You can use approxsoln (after correcting the B-problem pointed out by Carl):

AP:=NULL:
for k to 6 do
R := dsolve(eval({Eq1, Eq2, bcs1, bcs2}, B = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure,maxmesh=512,AP);
AP:=approxsoln=R;
X1 || k := rhs(R[3]); X2 || k := rhs(R[4]); Y1 || k := rhs(R[5]); Y2 || k := -rhs(R[6])
end do;

The print command at the end then gives me the values:

[-1.41421356296431, -1.43879570677334, -1.46322186668648, -1.48747814527757, -1.51155337283017, -1.55912751987572]

Incidentally, there is no reason for the 'print'. Just writing
[(X2 || (1 .. 6))(0)];

will do.