simplify does not work for sequences. See

simplify(sin(x)^2+cos(x)^2, 2*sin((1/6)*Pi)*cos((1/6)*Pi));

               Error, (in simplify/do) invalid simplification command

But

simplify([sin(x)^2+cos(x)^2, 2*sin((1/6)*Pi)*cos((1/6)*Pi)]);

                                        [1, (1/2)*sqrt(3)]

Replace the line   tab_values:=[evalf(simplify(seq(f(xx),xx=0..100)))];

by the line

tab_values := evalf([seq(f(xx), xx = 0 .. 100)]);

We prove by direct computation in Maple, that the point of intersection of any two perpendicular tangents to the ellipse  x^2/a^2+y^2/b^2=1   located  from its center at a distance  sqrt(a^2+b^2) . Suppose that the first coordinate of a point on the ellipse is  x1 .

restart;

f:=(x,y)->x^2/a^2+y^2/b^2-1:  F:=unapply(implicitdiff(f(x,y), y,x) ,x,y):

solve({f(x1,y1), f(x2,y2), F(x1,y1)*F(x2,y2)=-1}, {x2,y1,y2}):  # Finding the second coordinate of the first point and the coordinates of the second point with the perpendicular  tangent.

allvalues(%):

assign(%[1]):

solve({y=y1+F(x1,y1)*(x-x1), y=y2+F(x2,y2)*(x-x2)}, {x,y}):  # Finding the point of intersection of the two tangents

assign(%):

simplify(x^2+y^2);

                                                     a^2+b^2

x := xp-> (-1)*720.5668720*sinh(0.2043018094e-3*xp-0.8532729286)+84952.59423+4.014460003*10^5*arcsinh(0.1219272144e-1*sinh(0.2043018094e-3*xp-0.8532729286)-0.2032498888):

y:=yp->fsolve(x(xp)=yp, xp=0..10000):  # The inverse function

y(1),  y(100),  y(1000),  y(5000);  # Examples of use

                0.8684345847, 87.37094246, 923.0671034, 5469.479281

Edited.

For removal of parentheses in all the expression use  expand  command. For removal of parentheses in some parts of the expression, you can use  applyop  command.

Examples:

expand((a+b)^3);

applyop(expand, 1, (x+1)*(2*x-1)+(x-1)^3);

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 .