@Carl Love You forgot rhs at the end:

eval(intE, u= unapply(rhs(Sol),x));

@torabi 
Maybe I don't understand the problem, but to me the answers are simple:
1. The results are kept in the table res2. So pick the the one you want, e.g. res2[indx[1]].
2. Having found one omega2 you can try with a smaller or larger guess to find another.

Why don't you give us the Maple code you used? That may be easier to understand.

@sunflower p[-1] refers to the last element in p. When using sets you should be careful with referring to a particular element number, since it might not be what you thought. In this case it is no problem since Maple will sort integers in a set in increasing values.

Just try:
restart;
p:={13,2,9,5};
##You will see that Maple responds with
       p := {2, 5, 9, 13}

Compare this with the Maple response from executing the following definition of a list:

L:=[13,2,9,5];
##You will see that Maple responds with a copy in the same order
      L:=[13,2,9,5];

Your system can be seen as being the laplace transformed version of the following linear and homogeneous system of odes with initial values all equal to zero:

sys:=seq(diff(p[i](t),t)+mu[i]*p[4](t)=lambda[i]*p[0](t),i=1..4);
inttrans[laplace]~([sys],t,s);
ode:=diff(p[0](t),t)+add(lambda[i],i=1..4)*p[0](t)=add(mu[i]*p[i](t),i=1..4);
inttrans[laplace](ode,t,s);
SYS:=ode,sys; #The full system
dsolve({SYS,seq(p[i](0)=0,i=0..4)}); #Zero solution as expected
dsolve({SYS,seq(p[i](0)=0,i=0..4)},{seq(p[i](t),i=0..4)},method=laplace); #Ditto

@ Yes that was my point: Unless the 4 mu's are equal and the 4 lambda's are equal, there is only the trivial zero solution.

## Actually there is only the zero solution for any constants mu and lambda. No exceptions!

Your system is linear and homogeneous, so try some linear algebra:

A,z:=LinearAlgebra:-GenerateMatrix(eqs,Ps(s));
det:=LinearAlgebra:-Determinant(A);
##The determinant det has to be zero identically in s (at least for s>s0 for some s0).
solve(identity(det,s),{seq(mu[i],i=1..4),seq(lambda[i],i=1..4)}); NULL output
limit(det,s=infinity); # result 1

Your questions are not sufficiently clear to me to give you an answer.

The equations don't have any other solution than the trivial:
Ls := solve(eqs,{seq(P[i](s),i=0..4)});

unless all mu's are equal and all lambda's are equal. (Correction: no exceptions: see below.)

@Iza With the constants that you have given us so far:
mg:=2:l:=20:h:=5:xP:=l/5:P:=10:
l0:=40:

there is no root of eq1 in the interval 100..170.

So what has changed?

#####  But the problem is the order of the substitutions in your subs commands.
Notice that C has H in it:

has(C,H);
                              true
##Thus the FIRST substitution you make should not be H=H1. So use:
z1:=subs(C,H=H1,z1):
#z1:=subs(H=H1,C11=C[1],C21=C[3],z1); # WRONG order
#z2:=subs(H=H1,C12=C[2],C22=C[4],z2); #Ditto
z2:=subs(C,H=H1,z2):

## Since you are using Maple 16 you cannot use eval[recurse], which is easier to use because you don't have to think about the order:
z1:=eval[recurse](z1,C union {H=H1});

It is available in Maple2015 and Maple 18.

@want to be a permanent vegan For me the problem is the mathematics: If q stands for q(t) and the integration is really over q then the integral


must be rewritten as


where q(t1) = 0 and q(t2) = Pi.

If the integration is not over q, why on earth do you use q as variable of integration? Is it simply a mistake?
This is a mathematical question, not a Maple question.

In a reply above you write: "Because I don't know how to integral with the variable of t so I used q."
That I surely don't understand: You must first write down the integral you actually want whether you know how to integrate it or not.
Then you could ask for help on how to do it.

@want to be a permanent vegan 
In view of the familiar rule for changing variable:

A:=Int(f(q1),q1=a..b);   
          
IntegrationTools:-Change(A,q1=q(t),t);
   

you will need to change your integrals to include the derivative of q(t).
The limits of your q integrations I don't understand.

May we assume that the variable of integration in

sqrt(2)*(int(sqrt(H+w*cos(q)), q = 0 .. Pi, numeric))/Pi; #(Pi not pi)

is not supposed to be q, but t?
So that cos(q) would be cos(q(t))?

@Iza If C11, C12, .. depend on x, the integration cannot be done before substituting values, so my last suggestion certainly won't work.

It could be that your whole problem should be done numerically. To get any further you need to give us all the details.

### To illustrate what I have in mind here is a third version of the present problem.
## This version still assumes that C11, etc. are constants.
restart;
z1:=H/mg*cosh(mg/H*(x+C11))-C21;
z2:=H/mg*cosh(mg/H*(x+C12))-C22;
tg_alpha1:=diff(z1,x);
tg_alpha2:=diff(z2,x);
r1:=H*subs(x=xP,tg_alpha2)-H*subs(x=xP,tg_alpha1)-P;
r2:=subs(x=0,z1);
r3:=subs(x=l,z2)-h;
r4:=subs(x=xP,z1)-subs(x=xP,z2);
mg:=2:l:=20:h:=5:xP:=l/5:P:=10:
##We make a numerical procedure for computing eq1:
eq1:=proc(H1) local C,s1,s2,s,L0;
  C:=fsolve(eval({r1,r2,r3,r4},H=H1),{C11, C12, C21, C22}); #fsolve, H1 concrete
  s1:=evalf(eval(Int(sqrt(1+tg_alpha1^2),x=0..xP),{H=H1} union C)); #Numerical integration
  s2:=evalf(eval(Int(sqrt(1+tg_alpha2^2),x=l..xP),{H=H1} union C)); #Ditto
  s:=s1-s2;
  L0:=40;
  L0-s
end proc;

   
eq1(10); #Test
plot(eq1,10..11,adaptive=false,numpoints=10);
fsolve(eq1,10..11);

@Iza The integration is not the problem at all as you will see more clearly if you do the integrals before evaluating the constants C11 etc. for the values found by solve.
I have put that at the bottom of my answer.

Using "hello there" or `hello there` works for me. (I just tried: I never use spreadsheets).