@Carl Love You mentioned a named expression sequence and gave an example of a list of lists.
So consider:
f:=(A, index::list)-> A[index[]];
LoL:=[[2,3],[8,9]];
map(f,LoL,[-1]); #Works as intended
map(`?[]`, LoL, [-1]); # So does this
SoL:=[2,3],[8,9]; #Named sequence of lists
map(`?[]`, SoL, [-1]); #As expected it doesn't work: Error
map(f, SoL, [-1]); # Doesn't work as maybe intended
## But the following works for `?[]`, but not for f:
`?[]`(SoL,[-1]); #Works
f(SoL,[-1]); #Error

To me this sounds like the familiar shooting method.
Do I understand you correctly if I interpret your statement
   "should work the package of optimization in relation to a system of nonlinear equations"
as meaning that instead of using fsolve (or solve) on the equation y(1) = Y1 you would minimize (e.g.) the (sum of) square(s) y(1) -Y1 ?
That would just amount to using another solver than fsolve, e.g. Optimization:-LSSolve.

@snhaseela2 Solving boundary problems for ode systems can be more challenging than solving initial value problems. So any method that works is worth considering. 

When successful dsolve/numeric/bvp is easier to use.

With Douglas Meade's Shoot package you need to rewrite the system as a first order system and you have to provide a start guess. Finding one that works may not be easy.
With the shooting method I gave you don't have to rewrite the system, but for success you need to supply a start guess to the solver used, in my case fsolve. The method as presented in my answer is not written as a ready made package and may be difficult to understand for a new Maple user. Writing a package is not difficult though. I did, but I'm not convinced that the result is good enough for public consumption.

Even dsolve/numeric/bvp in fact needs an initial "profile". It tries to find one itself. If unsuccessful the user can supply an approximate solution, approxsoln=[f(eta)= ...., theta(eta)= ..., phi(eta)= ... ].

@snhaseela2 You have several problems of syntax:

1. All unknown functions (A,B,C, etc) must appear as A(eta), B(eta), C(eta) etc.

2. You need to use another name than D (e.g. D1) since D is used as the differention operator on functions: Thus D(sin) = cos.
3. You cannot use square brackets, [], instead of parentheses, (), for grouping.
4. You cannot leave out some multiplication sign. In 2D math it is enough with a space, but my advice is to use * always.
As an example you have MB without space or *.

Also clearly beta needs to be given a value.
All these syntax problems must be dealt with. Take a close look at eval(ODE, beta = 1) before proceeding: Does it look right?

Probably what you are using is the shoot procedure made by Douglas Meade, which (reasonably enough) needs a start guess. Try the one I used in my answer - and try my answer too!

@snhaseela2 You may not be aware of odeplot from the plots package. It is very useful.
Here is how it could be used in your case using Tom Leslie's revised version of yours.
Notice that I keep all the 4 results in their entirety (as res[1], res[2] etc. )

for k to 4 do
  res[k]:= dsolve(eval({Eq1, Eq2, Eq3, bcs1}, beta = L[k]),numeric, method = bvp[midrich], maxmesh = 4096)
end do;
#Now plot whatever function you like:
plots:-odeplot(res[1],[[eta,D(f)(eta)],[eta,theta(eta)],[eta,(D@@3)(f)(eta)],[eta,-phi(eta)]],0..5);
#Here is a short animation in beta:
S:=seq(plots:-odeplot(res[k],[[eta,D(f)(eta)],[eta,theta(eta)],[eta,(D@@3)(f)(eta)],[eta,-phi(eta)]],0..5),k=1..4):
plots:-display(S,insequence);

##For information about maxmesh go to ?dsolve,numeric,bvp
##The method used by dsolve is not a shooting method. A shooting method formulates the problem as an initial value problem at some point with some given values and some variable ones. Then it tries to find the correct values of the variable ones. In your case you could try shooting from eta=blt (not eta=0 since f(0)=0 is imposed and f(eta) appears in the denominator for the fourth derivative of f. Try
solve({Eq1, Eq2, Eq3},{diff(f(eta),eta$4),diff(theta(eta),eta,eta),diff(phi(eta),eta,eta)});

You define N twice, but don't give a definition of M.
By the way: Use Vector instead of the deprecated vector.

@iman You are plotting w.r.t. sigma1, but sigma1 is set at 40. Did you mean sigma2?

I didn't try solve for very long, so never got a result from that. Switching to fsolve I tried

S:=proc(s2) local res;
   res:=fsolve(eval({Q1, Q2, Q3, Q4},sigma2=s2));
   if not type(res,set) then print(cat("failure at sigma2=",s2)) else eval(p1^2+q1^2,res) end if
end proc;
Digits:=15:
plot(S,-100..200,adaptive=false,numpoints=100);

@Markiyan Hirnyk The title of your comment was:  No copyright.

@Markiyan Hirnyk What is the point of saying that the print from

print(`evala/toprof`);

shows nothing but proc(a) ... end proc ?

As you well know by using

interface(verboseproc=2);

you get the whole deal.

And did you try the command given by Dave L:
op(3,eval(`evala/toprof`));

@iman Use eval (as here) or subs:

res := solve({pprime11, qprime11}, {p1, q1});
nops([res]); #Answer 2
res[1]; #Uninteresting I suppose, so we use res[2]
S:=eval(p1^2+q1^2,res[2]):
Digits:=15:
plot(S,sigma1=-50..50);

@iman OK, I get your point: the unknowns are p1 and q1.

As pointed out by "one man" the solution involves RootOf. That is for the simple reason that you need to find the roots of a polynomial of degree 8:

res:=solve({pprime11,qprime11 },{p1,q1});
pol:=op([1,1],indets([res],specfunc(RootOf))); #The polynomial (in _Z)
degree(%,_Z); # Result 8

Your two expressions are cubic polynomials in sigma1:
degree(pprime11,sigma1); #result 3
degree(qprime11,sigma1); # result 3

The 3 roots of each are found by
res_p := solve(pprime11, sigma1);
res_q:=solve(qprime11,sigma1);

What is your intention? To determine p1 and q1 such that this pair of 3 roots are identical?

@Doug Meade I get
D(f)(3/8) = 1/2+(1/4)*sqrt(10) (or D(f)(3/8) = 1/2-(1/4)*sqrt(10) )

by doing
eval[recurse](convert(deq1,D),{b=3/8,f(3/8)=0});
solve(%,{D(f)(3/8)});


Another comment: f''(b) vanishes from deq1 at b = 3/8 under the assumption that it stays bounded in a neighborhood of b = 3/8.
A priori it is possible that b*f''(b) has a nonzero limit as b -> 3/8 from the right.

Just a couple of observations:

Didn't you mean delimiter = " " , i.e. with a space inside "" ?

I noticed that without datatype = string the first works:
ImportMatrix(LambdaFile, delimiter = " ");

##Furthermore, if the line change to an empty line is removed from the first file then also this works:
ImportMatrix(LambdaFile, delimiter = " ",datatype=string);

Your second file also has a line change to an empty line, but as you say there is no problem.

@Carl Love I didn't try with matrices. But now I get your point.
Actually || beats cat in this case:

restart;
for i from 0 to 2 do cat(CCnew,i):=LinearAlgebra:-RandomMatrix(3) end do;
`<|>`(-CCnew||(0..2)); #Doesn't work as intended
## A two step procedure works:
-CCnew||(0..2);
`<|>`(%);
## So this should work:
`<|>`(eval(-CCnew||(0..2))); #Works
## Now cat:
`<|>`(eval(cat(-CCnew,0..2))); #Doesn't work as intended at all (with or without eval)
lprint(%[1]); # !!! symbol :  `-CCnew0` # Since Maple 2015. cat was updated then.
## For cat we need eval even when there is no minus:
`<|>`(cat(CCnew,0..2));
%; ##Two step OK
`<|>`(eval(cat(CCnew,0..2))); #Works