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.

@torabi A few comments:

Since you say that some result is not correct you must know what "the correct result" is. How do you determine what result is correct?
As an example, what would be the correct value of the parameter a in my pendulum example above? My answer would be that they are all correct. a is closely related to the period of the pendulum which depends on the amplitude in the nonlinear case.

Your worksheet "first result .." has a linear system (in effect) as I already mentioned: When s(x) = 0 for all x then the equation for g is linear.

@torabi I'm convinced that you are wrong. The values for omega2 are really different.

Consider the familiar ode for the pendulum for a simpler example:

restart;
ode:=diff(y(x),x,x)+a*sin(y(x))=0;
bcs:={y(0)=0,y(1)=0};
res:=y1->dsolve({ode} union bcs union {D(y)(0)=y1},numeric);
res(5)(0); # a = 13.3616633336330
res(0.1)(0); # a = 9.87085448441260
#Try the lines above with
#ode:=diff(y(x),x,x)+a*y(x)=0;
# and you will find the same values of a in the two cases (except for roundoff errors).
## Just for animation purposes I define AY:
AY:=proc(y1) option remember; local res; 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);
  eval(a,res(0)),res
end proc;
plots:-animate(plots:-odeplot,[AY(y1)[2],[x,y(x)],0..1,caption=typeset(a = AY(y1)[1]) ],y1=0.1..5);

@torabi My point is that to each value (s1) of s(1) corresponds a value of omega2. Thus omega2 is a function of s1. In the linear case that function is constant, not so in the nonlinear case. If your intention is to find several (all) values of omega2 you therefore have a problem. What you can do is to find pairs [omega2, [s,g] ] consisting of omega2 and a corresponding solution [s,g]. On the surface that sound like the situation in the linear case, but in that case to omega2 corresponds not only a particular [s,g], but also all constant multiples.
So a rhetorical question to you could be: Which of the omega2 values found by my code above with
s(1)=-.3+0.051*i , i = 0 .. 12, are you going to accept as the one you want?

@torabi It may be of interest to try other values for s(1) when using the s(1) = something as the extra boundary condition. For a linear system this wouldn't be interesting since all constant multiples of a solution (s,g) corresponding to some omega2 would again be a solution with the same omega2. But for a nonlinear system this is no more the case.
So I tried the following, where I used try...catch..end try despite my comment above.
Define sys and bcs as in my answer above.
Digits:=15:
AP:=NULL:
for i from 0 to 12 do
  try
    #To avoid s(1)=0 we choose the increment 0.051 rather than 0.05
    res[i]:=dsolve(sys union bcs union {s(1)=-.3+0.051*i},numeric,AP,abserr=1e-10);
    AP:=approxsoln=res[i];
  catch:
    AP:=NULL;
  end try
 end do;
ind:=sort([indices(res,nolist)]);
plots:-display(seq(plots:-odeplot(res[i],[[x,g(x)],[x,s(x)]],0..1),i=ind));
plots:-display(seq(plots:-odeplot(res[i],[[x,g(x)],[x,s(x)]],0..1,thickness=3),i=ind),insequence);
##The values for omega2:
seq(subs(res[i](0),omega2),i=ind);
sqrt~([%]); #The values for omega
sort(%);
evalf[5](%);
#Question: Are these values of omega really different? I think so, and why not, the system is nonlinear.

@torabi Let me answer your questions 3 and 4 first:

odetest({s(x) = 0}, sys)  is done to see under which circumstances sys has a solution with s(x) = 0 for all x. This turns out to be the requirement that g satisfies the linear ode:
 -0.326905829596411e-2*g(x)-(diff(g(x), x, x))-(4/3)*omega2*g(x) = 0.
So the answer to question 4 is to look at the value of the coefficient of g(x) in that ode (I called it ode_g).
But notice that 3 and 4 deals with the very special case where s(x) is identically zero. This case may not be of physical interest (assuming that this is an applied problem). For that case all eigenvalues have been determined.

For the general case where s(x)<>0 you have to try various nonhomogeneous extra boundary conditions one at a time just as I did with the example
res:=dsolve(sys union bcs union {s(1)=.1},numeric);
Just replace s(1) = .1 with your choice. You may not succeed with that choice, then try another.

I once gave an answer to a similar problem in MaplePrimes where I used a loop with a try ... catch .. end try structure looping through simple choices of an extra boundary condition. The try ... catch .. end try structure would catch any error and if so the loop would continue with the next choice. I'm not going to repeat that here since it turned out that it was harder to understand and therefore more difficult to adapt to other circumstances.

To find other values of omega2 (if there are any) try using an approximate solution with more wiggles.
approxsoln=[s(x)= something, g(x) = something]. See the help page for dsolve,numeric,bvp.

@Markiyan Hirnyk Quite correct. Notice the closing remark in my answer.

@Carl Love I got the same error.
I just tried to isolate the highest derivatives:
solve({Eq1, Eq2, Eq3, Eq4},diff~({f,theta,h,g}(eta),eta,eta));
## That suggested using a midpoint method (  sqrt(diff(f(eta), eta)) being in the denominator):
dsolve(eval(dsys,lambda=0), numeric,method=bvp[midrich]);
## which didn't solve the problem.
## The imaginary I may come from the sqrt in  sqrt(diff(f(eta), eta)) (I changed nn from 1.5 to 3/2).
## I shall have another look in a few hours.

The multiplication operator in Maple is * not . (dot). So begin by changing all those dots to *.
After that you still have a problem.
I suggest omitting continuation in lambda until you have a problem that only suffers from convergence issues.

@mskalsi If you are using Windows try unzipping the mla file to:

"C:/Users/xxxx/maple/toolbox/DifferentialGeometry"

You need to replace xxxx with the name of the user (presumably you). Don't change anything else!

Then in the worksheet supplied by Ian Anderson (Optima1.mw) start with

restart;
libname := "C:/Users/xxxx/maple/toolbox/DifferentialGeometry", libname;
### Then continue with
with(DifferentialGeometry);
## as in the worksheet.
## Using that location should ensure that the help also works: Try
?DifferentialGeometry
## I just downloaded and tried the above with xxxx=Bruger (which happens to be my username on this computer). It worked.