Maybe OP meant the creation of new functions from existing functions by using the composition of the operators and other operators?

An example:

x:=t->t^2:   y:=t->sin(t):

f:=x@(y+x):

g:=y@(x-y):

f(t);

g(t);

To find the values of a function better  to specify it as a function not as an expression:

g:=y->200*y:

g(4);

              800

The workaround for your method:

g(0.5,y):=200*y:

eval(g(0.5,y), y=4);

               800

From the standpoint of rigorous mathematics, evalf  command is not sufficiently reliable. Here is a simple example of a divergent series and finding its "sum" with  evalf :

f:=unapply(`if`(n<10^10, 1/(2^n+3^n), 10^(-10)), n);

evalf(sum(f(n), n=1..infinity));

Obviously, Maple ignores the behavior of series for  n>=10^10 , and calculates the sum by using the first few terms.

For a rigorous solution of the problem can make the following two steps:

1) The proof of the convergence of the series easy to make, using  d'Alembert's ratio test .

2) For a rigorous calculation with a given accuracy is necessary to make an estimate of n-th remainder of the series R[n] = S - S[n] (here  S  is the sum of a series and  S[n]  is  n-th partial sum).

Step 1 we make with Maple:

u:=n->(n!/(2*n)!)^n:

limit(u(n+1)/u(n), n=infinity);

                                0

Step 2 by hand:

First we estimate  u[n+1]  and then estimate  R[n]  by an geometric series  sum(b*q^n, n=0..infinity)=b/(1-q)  with  q=1/2

Obviously,   

u(n+1)=((n+1)!/(2*(n+1))!)^(n+1)=(1/((n+2)*...*(2*n+2)))^(n+1)

<(1/(n+2))^(n^2+2*n+1)

This is the final estimate:   R[n] < 2/((n+2)^(n^2+2*n+1))

Example of use.

It is necessary to calculate the sum of the original series with an accuracy of 10^(-10) :

fsolve(2/((n+2)^(n^2+2*n+1))=10^(-10));

                        2.870597417

We see that in order to achieve the given accuracy is sufficient to find sum the first three terms of the series.

See the  result

In the assignment  

f:= x(t)^2 + y(t)^2:  

f  is an expression, not a function. This is inconvenient, if necessary to calculate the value of the function and its derivative at a point. See

f:= x(t)^2 + y(t)^2:

diff(f, t):

f(1);

                                   x(t)(1)^2+y(t)(1)^2

diff(f, t)(1);

                   2*x(t)(1)*diff(x(t),t)(1)+2*y(t)(1)*diff(y(t),t)(1)

It is better to define  f  and its derivative as  functions:

f:=t->x(t)^2+y(t)^2:

D(f);

                        t -> 2*x(t)*D(x)(t)+2*y(t)*D(y)(t)

f(1);

                                      x(1)^2+y(1)^2

D(f)(1);

                         2*x(1)*D(x)(1)+2*y(1)*D(y)(1)

Specific example:

f:=t^2;  # Assignment to the expression

                                      f := t^2 

f(1);  # Here Maple does not understand that it is necessary to substitute the number  1  instead of  t

                                       t(1)^2

diff(f, t);

                                         2*t

diff(f, t)(1);  # The same for the derivative

                                        2*t(1)

f:=t->t^2;  # The functional assignment

                                    f := t -> t^2

D(f);  # The derivative as a function

                                        t -> 2*t

f(1);  # Maple correctly calculates the value of the function  f  at the point  t=1

                                           1

D(f)(1);  # The same for the derivative

                                           2

 Edited.

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