As a practice in memorizing the alphabet I tried the following variant of the first example (with unevaluation quotes as in my previous comment):

restart;
expr:='a+b+c+d+e+f+g+h+i+j+k+l+m+n+o+p+q+r+s+t+u+v+w+x+y-y+z';
#sort(expr); #See comment below
expr; #Correct
expr; #Not correct and why are the two z's not combined as 2*z ?
simplify(%); #No effect
indets(%,name); #The whole alphabet
nops(%); # 26
### If you run the code with the sort command uncommented, then the result of sort(expr) is correct, but the rest is not.
## Try continuing with:
eval(expr,{z=Z,y=Y}); #Result: the two separate z's have been replaced by 2*Z and y by Y.


The same can be done with your symbols:
restart;
expr:='t1+t2+t3+t4+t5+t6+t7+t8+t9+t10+t11+t12+t13+t14+t15+t16+t17+t18+t19+t20+t21+t22+t0-t0+t23';
expr; #Correct
expr; # Incorrect

I tried putting unevaluation quotes around expr1 in your very first example:

expr1:='t1+t2+t3+t4+t5+t6+t7+t8+t9+t10+t11+t12+t13+t14+t15+t16+t17+t18+t19+t20+t21+t22+t0-t0+t23';

The response was the very same expression without the quotes. Thus automatic simplification of t0-t0 didn't happen. That I find strange too.
With expr2 defined exactly the same:

expr2:='t1+t2+t3+t4+t5+t6+t7+t8+t9+t10+t11+t12+t13+t14+t15+t16+t17+t18+t19+t20+t21+t22+t0-t0+t23';

you get exactly the same as for expr1, but we still have

expr1-expr2;
                         t0-t23

The addresses of expr1 and expr2 (in this case with unevalution quotes) are the same:
addressof(expr1) - addressof(expr2); #zero

@want to be a permanent vegan The response from my browser when clicking your link  hw2....mw is "Bad Request".

But I must insist that it is a mathematical question: Formulate your problem mathematically, then there are quite a few people in this forum, who could help you with the code.

I don't know why you don't take my previous comments to a related question seriously.

http://www.mapleprimes.com/questions/206784-How-To-Achieve-A-Piecewise-Function

In your new worksheet you have the integral

int(sqrt(2*(H(t)+omega*cos(q(t)))), q(t) = q(t)-2*Pi .. q(t));

which presents the same problems I mentioned there. I even gave the solution.

I shall take a look at it.
In the meantime I think you should replace Q, Q1, Q2, Q3, Q4 by the following. They do the present job simpler (and faster) than the one you got (for which I am to blame):

Q:=proc(pp2,fi0,HTCC,ti0) option remember; local res,F0,F1,F2,F3,a,INT0,INT10,INT15,INT20;
print(pp2,fi0,HTCC,ti0);
if not type([pp2,fi0,HTCC,ti0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({subs(phi0=fi0,p2=pp2,HTC=HTCC,eq2)=0,subs(phi0=fi0,eq1=0),eq3=0,D(u)(1)=0,u(0)=lambda*D(u)(0),(T)(1)=ti0,phi(0)=fi0,D(T)(1)=0}, numeric,method=bvp[midrich],output=listprocedure);
F0,F1,F2,F3:=op(subs(res,[u(eta),phi(eta),T(eta),diff(T(eta),eta)])):
INT0:=evalf(Int(2*F0(eta)*(1-eta),eta=0..1));
INT10:=evalf(Int(2*F0(eta)*F1(eta)*(1-eta),eta=0..1));
INT15:=evalf(Int((2*F0(eta)*(1-eta)*(F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )),eta=0..1));
INT20:=evalf(Int((2*F0(eta)*F2(eta)*(1-eta)*(F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )),eta=0..1));
a[1]:=INT15-10000*pp2;
a[2]:=INT10/INT0-phi_avg;
a[4]:=F2(0)-1:
a[3]:=INT20/INT15:
[a[1],a[2],a[3],a[4]]
end proc;
Q1:=proc(pp2,fi0,HTCC,ti0) Q(_passed)[1] end proc;
Q2:=proc(pp2,fi0,HTCC,ti0) Q(_passed)[2] end proc;
Q3:=proc(pp2,fi0,HTCC,ti0) Q(_passed)[3] end proc;
Q4:=proc(pp2,fi0,HTCC,ti0) Q(_passed)[4] end proc;

@Markiyan Hirnyk Yes, I tried something similar. I used
dsolve({Eq,D(y)(0)=0});
with the intention of getting a solution which is an even function. Indeed you get an even function it appears. The value at x=0 is 0, but the value at x=1 and at x=-1 is (.444221192534724+0.*I)*K (when using evalf on the result).

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