Unfortunately, plots:-inequal  command for nonlinear inequalities and   plots:-shadebetween  command work only in the latest versions of Maple. For older versions, you can use  plots[implicitplot]  command with  filledregions  option .

For example in Maple 12:

plots[implicitplot]((y-sin(x))*(y-cos(x)), x = 0 .. 4*Pi, y = -1 .. 1, coloring = [red, white], filledregions = true, gridrefine = 5, scaling = constrained);

 In new versions this method also works.

which is built on the increase of the parameter  T :

plots[animate](plot, [[cos(t), sin(t), t = 0 .. T], color = red, thickness = 2], T = 0 .. 2*Pi, scaling = constrained, view = [-1.2 .. 1.2, -1.2 .. 1.2], frames = 60);

The upper axis is made artificially by the objects  B, C, T :

restart;

A := plots[dualaxisplot](x^2, x^2, x = 0 .. 4):

B := plot(16, x = 0 .. 4, color = black, thickness = 0):

C := seq(plot([i, t, t = 16-0.15 .. 16], color = black, thickness = 0), i = 1 .. 3):

T := plots[textplot]([seq([i, 15.5, i], i = 1 .. 3)]):

plots[display](A, B, C, T);

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);