@acer You are right (of course).
I just checked Fede's worksheet without making any changes and using Maple 18.02 on my old Windows 7, 32 bit machine. The timings were roughly the same. I encountered no problems at all with NumSteps= 6, 20, 40, 80.

It would help if you uploaded a worksheet. Use the green arrow in the MaplePrimes editor.

@acer Using method = _d01amc sure speeds up the computation: In Maple 18.02 (with the matrix<->Matrix replacement I mentioned) I got 0.44 seconds on my old machine and 0.17 on my new machine with Maple2015.2.

You use matrix instead of Matrix in the definition of A. That may make a difference. (Apparently not: see below).
I changed to:

A := Matrix(NumSteps, 2, datatype = float);

Similarly, I added the datatype option to the definitions of T1 and AF.

I measured the timing of the loop only.
For NumSteps=40 in Maple 18.02 on an old not very fast machine I got 21 seconds. For NumSteps=80 the time used was 44 seconds.

In Maple 2015.2 on a new and reasonably good machine the times were about half these amounts.

@want to be a permanent vegan You refer to 3 worksheets. I had a brief look at one of them (example_hw2.mw) , which has these integrals:

and

Since p(t) apparently doesn't depend on x these two integrals are obviously easily computed, just as in

restart;
int(p(t),x=-eta(t)..eta(t));
                             

Surely you don't intend anything as trivial as that?
So I'll stick to my statement about mathematics: Formulate the problem mathematically, so it is obvious to any person with some background in mathematics.
The point is (as somebody said): If you don't know where you want to go, any road will take you there.

@crazyeti You wrote

"I just hope Maple could post a solution/workaround soon, something like always executing the assignment statement twice etc."

That won't work as it is the second one that is wrong.

Added: I have no idea about what is going on, but it is somewhat entertaining (Maple2015.2):

restart;
expr:='a+b+c+d+e+f+g+h+i+j+k+l+m+y-y+n+o+p+q+r+s+t+u+v+w+x+z';
expr;
expr;
indets(expr,name);
nops(%);
nops(expr);
eval(expr,{w=W,z=Z}); #2*W and 2*Z

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.