My first reaction is that a thing like

(D@@2)(phi)(theta(t))(t); #(I removed superfluous parentheses)

is somewhat weird.

It means the second derivative of phi taken at theta(t) and then evaluated at t. That must be rethought!

1. You call the expression exp(-h^2/T)*(Int(exp(-x^2/T)*BesselI(0, h*x/T)*x, x = 0 .. 1))/T an equation. So is that expression supposed to be equal to zero?

2. What is the unknown: h or T?

4. BesselI(0,y) appears to be positive for all real y, so it seems you have a problem.

So what do you expect us to do?

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