@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

@Carl Love By right next to CCnew do you mean ((-CCnew)||(0..dmax-1) ? (i.e. with parentheses around -CCnew.)
I tried the following (where the use of `<|>` is actually irrelevant):
restart;
C:=A;
`<|>`(-C||(0..6)); #No problem
`<|>`((C)||(0..6)); # Problem: Error, `||` unexpected
## || is expecting a name or a string.
## The first argument is not evaluated by ||, but the problem probably is at the parsing level
## The help page for || discourages the use of ||. cat can be used instead:
`<|>`(cat(-C,0..6));
`<|>`(cat((C),0..6));
`<|>`(cat(-'C',0..6));

@torabi In my comment "Correct result?" I wrote

"a is closely related to the period of the pendulum which depends on the amplitude in the nonlinear case."

I was referring to the pendulum. I asked you which of the values found for a you would say is correct. You haven't told me. In fact of course they are all correct: The period of the pendulum is dependent on the amplitude.
But for small amplitudes it is roughly constant. For small amplitudes (angles y) we can replace sin(y) with y.
For the pendulum I found only the lowest values of a corresponding to solutions with no zero in the interval.
To get the solutions with one zero we can use an approximate solution with one zero:

res2:=y1->dsolve({ode} union bcs union {D(y)(0)=y1},numeric,approxsoln=[a = 40, y(x)=sin(2*Pi*x)] );
res2(5)(0);
res2(0.1)(0);
AY2:=proc(y1) option remember; local res,aa; global ode,bcs,a;
  if not type(y1,realcons) then return 'procname(_passed)' end if;
  res:=dsolve({ode} union bcs union {D(y)(0)=y1},numeric,approxsoln=[a = 40, y(x)=sin(2*Pi*x)]);
  eval(a,res(0)),res
end proc;
plots:-animate(plots:-odeplot,[AY2(y1)[2],[x,y(x)],0..1,caption=typeset(a = AY2(y1)[1]) ],y1=0.1..5);

Again we notice the dependency of a on the amplitude characteristic of the nonlinear case.

Now for your example you use b = 0.001 in the case Omega=30 and b = 0.01 when Omega = 0, where b runs through
extra_bcs = {g(1), s(1), (D(g))(0), (D(g))(1), (D(s))(1), ((D@@2)(s))(0), ((D@@3)(s))(0)}.

Maybe for your system the values found for omega2 are roughly constant for the small values of the successful b's in extra_bc, but you have to look at the graphs also as there may be other ranges of omega2 corresponding to more zeros as I illustrated with the pendulum. 

I took a brief glance at the pdf-file you uploaded, but cannot spend any time on that, sorry.