You can do any branching conditions, if you are using  a nested  `if`.

Example:

expr:=`if`(x>-3 and x<-2,6*(x+3),`if`(x>=-2 and x<0,-3/2*x+3,`if`(x>=0 and x<2,3/2*x+3,`if`(x>=-2,-6*(x-3),undefined)))):

 plot(expr, x=-3..3, color=red, thickness=3);

b:=[3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71];  # b is defined as a list

exp(1) -~ b;  # elementwise subtraction

evalf(%);

I converted your code into 1D-math. At first

beta__1 := 2*U__n/U;

...

but then

Tmaxinverse := int(2*beta[1]*(-y^2+1)*(1+beta[2]* ...

beta__1  is not equal to  beta[1]

 Replace  beta__1 := 2*U__n/U;  with  beta[1] := 2*U[n]/U;   and so on

It is even easier:

restart;

f:=(x,y)->evalf[5](2*x*Int(sqrt(1+y^2*(t*x-1)^2/(1-(t*x-1)^2)), t = 0 .. 1)):

printf(" x     y     f(x,y)\n");    

for x from 0.1 to 1.9 by 0.1 do

for y from 0.1 to 0.9 by 0.1 do

printf("%g   %g   %g\n",x,y,f(x,y));

od; od;

You can output these values in the form of a matrix as follows: 

restart;

Digits:=5;

f:=(x,y)->2*x*evalf(Int(sqrt(1+y^2*(t*x-1)^2/(1-(t*x-1)^2)), t = 0 .. 1)):

for i to 20 do

for j to 11 do

a[i,j]:=f(0.1*i, 0.1*(j-1));

od: od:

interface(rtablesize=infinity):

A:=Matrix(20,11, (i,j)->a[i,j]):

V:=<seq(i,i=0.1..2,0.1)>:

W:=Vector[row]([cat(x,` \\ `,y),seq(i,i=0..1,0.1)]):

<W, <V|A>>;

1) In the syntax (in your file)  both ways of calculating the integral are wrong because you wrote that function  p[t]  depends on  theta  and  r .  In calculating the integral Maple believes that  p[t]  does not depend on  theta  and  r . 

2) If you write explicitly the function  p[t](theta,r)  and  r[5], r[6]  and  theta[tt]  are some constants, it is possible to use both methods. Their efficiency is apparently dependent on the form of the function  p[t](theta,r).

You have an inequality with 3 variables. In general, the solutions will be some spatial regions. Maple  has not any simple description of these regions. But it is easy to find symbolic (or numerical) solution if the two variables are specified.

Example:

Ineq := 3*a[1]*a[2]*a[3]-2*a[1]*a[2]-2*a[1]*a[3]-2*a[2]*a[3]+2*a[1]+2*a[2]+2*a[3]-1 <= 3*a[1]*a[2]*a[3]-((a[1]*a[2]*a[3]+1)/(a[2]*a[3]-a[3]+1)-1)*((a[1]*a[2]*a[3]+1)/(a[1]*a[3]-a[1]+1)-1)-((a[1]*a[2]*a[3]+1)/(a[2]*a[3]-a[3]+1)-1)*((a[1]*a[2]*a[3]+1)/(a[1]*a[2]-a[2]+1)-1)-((a[1]*a[2]*a[3]+1)/(a[1]*a[3]-a[1]+1)-1)*((a[1]*a[2]*a[3]+1)/(a[1]*a[2]-a[2]+1)-1):

solve(eval(Ineq, {a[1] = 3, a[2] = 5}));

simplify(fnormal([evalf(%)]), zero);

eq1 := -3.999168585*10^5*p1*q1^2-2.543619525*10^7*q1*p1^2+7.762255614*10^6*q1^3-3.999168584*10^5*p1^3-0.5000000000e-4*p1+0.1007754033e-3*q1 = 0:

eq2 := -3.999168585*10^5*q1*p1^2-7.762255614*10^6*p1^3+2.543619525*10^7*p1*q1^2-0.5000000000e-4*q1-3.999168584*10^5*q1^3-0.1007754033e-3*p1 = 0:

solve({eq1, eq2});

I did not understand that  pseudo rule means .

bcount:=n->nops(select(t->type(t, odd), [seq(binomial(n,k), k=0..n)])):

Example:

bcount(6);

                             4

It would be interesting to find an explicit formula for this.

Edited.

For any simple polygon (not necessarily convex) area can be found by shoelace formula . The algorithm implemented in Maple in the procedure  Area, which also solves  more general problem - symbolically (or numerically) finds the area of a figure, bounded by a non-selfintersecting piecewise smooth curve. 

convert(sin(x), FormalPowerSeries);

Should be

restart;

Ec := (1/2)*Kc*(2*Pi*h0/Pi/(a0^2+h0^2)-2*Pi*h/Pi/(a^2+h^2))^2;

simplify(Ec);

You have a system of 4 matrix equations in 4 unknowns matrices and matrix parameter y . I think that Maple can not solve the system symbolically, but Maple can solve it numerically for specific values of  B2 , B3  and  y (I took all the square matrix of the third order) . 

I have replaced the names of all the matrices with capital letters and converted the original system in 36 nonlinear scalar equations.

restart;

with(LinearAlgebra):

local D;

B2, B3, Y:=seq(RandomMatrix(3, generator=0..10), i=1..3);

A,B,C,D:=seq(Matrix(3, symbol=i), i=[a,b,c,d]); 

eq2 := A.B+C.D+A:

eq3 := A.C+C.D+C:

eq4 := A.B+C.A+B.C:

eq5 := A.B+A.D+B.C:

fsolve(map(t->Equate(op(t))[], [eq2=B2,eq3=B3,eq4=B2,eq5=Y]));

assign(%):

A,B,C,D;

Expr:=expand((x^2-x-3)^10);

select(t->sign(t)=-1, [op(Expr)]);

nops(%);

restart;

solve({-(1/3)*(eta+5)/(eta-3) = 3*eta/(-2+eta)});

assign(evalf(%[1]));

eta;