I agree with Kitonum about solving numerically.
But I'd like to make a couple of observations:

Notice that if y(x) is a solution of Eq then so is y(-x). This doesn't imply that all solutions are even functions of x, but if your boundary value problem has a unique solution then that solution would be even.

Eq := diff(y(x), x, x) = -(x^2+1)*y(x)+K;
PDEtools:-dchange({y(x)=w(xi),x=-xi},Eq,[w,xi]);
#Same ode as output
##We try imposing just one condition:
res:=dsolve({Eq,y(1)=0});
##Notice the arbitrary constant _C2 in res.
##Let us procede with an attempt to determine _C2 from imposing y(-1)=0
y1:=eval(rhs(res),x=-1);
indets(y1,name);
#Notice that _C2 is gone! Thus there is no arbitrary constant we can use to make y(-1)=0.
#We should try if y1 is actually zero:
evalf(y1);
#You will notice that it is not.
The solution res also has an immediate problem at x=0:
eval(rhs(res),x=0);
##Taking limits instead:
LR:=limit(rhs(res),x=0,right);
LL:=limit(rhs(res),x=0,left);
##Those limits are not the same, but can be made so by choosing _C2 from that requirement.
solve(LR=LL,{_C2});res2:=eval(res,%);
#So now is y(-1) = 0 although we have no arbitrary constant to use to make it so?
evalf(eval(rhs(res2),{K=1,x=-1}));
#No!

@hind Maple has no facility to change the order of integration and surely for a good reason.

It is simple for a case like the one we have here, where the double integral is over a region, which can be described as the region between the graphs of two (simple) functions of x and also between the graphs of two equally simple functions of t.
If we remove the 'simple' in the statement above the change can still be made, but given the two functions of x, how to find the two functions of t in a programmatic way (assuming that they exist)?
Consider for an only slightly more complicated example the two bounded regions bounded by y=x and y=4*x*(1-x) on this picture:

plot([x,4*x*(1-x)],x=0..1,y=0..1);

The region bounded above by the parabola and below by the line y = x is indeed the region between the graphs of simple factions of x. Turn your head 90 degrees clockwise and you will see that if you replace "x" by "t" and remove "simple" you can say the same. The turned head upper graph is the graph of a piecewise defined function, not a major problem, but still a complication.
For the region above the x-axis and below both the line and the parabola the situation is reversed: The upper curve is the graph of a piecewise defined function. The turned head upper or lower are graphs of simple functions.

@Carl Love In this case all I'm saying is that

Int(Int(f(x,t),t=0..x),x=0..y)=Int(Int(f(x,t),x=t..y),t=0..y);
 
Each of the iterated integrals computes the double integral over a triangle T with vertices (0,0), (y,0), (y,y) given by
Int(f(x,t),A=T..``);



The triangle can be described as the domain above the x-axis and below the curve t=x, where 0 <= x <= y.
If you turn your head 90 degrees clockwise you see that it can also be described as the domain below the line x=y and above the curve x=t, where 0 <= t <= y.

@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))?