T:=table([(25, 1) = -39, (16, 151) = 32, (33, 1) = -54, (1, 1) = 29, (13, 1) = 32, (31, 101) = -7, (6, 51) = -10, (11, 101) = -1, (28, 151) = -39, (18, 51) = -65, (4, 151) = 29, (8, 151) = -10, (23, 101) = 23, (34, 51) = -54, (40, 151) = 87, (36, 151) = -54, (9, 1) = -1, (37, 1) = 87, (21, 1) = 23, (14, 51) = 32, (22, 51) = 23, (20, 151) = -65, (27, 101) = -39, (3, 101) = 29, (19, 101) = -65, (24, 151) = 23, (32, 151) = -7, (30, 51) = -7, (38, 51) = 87, (7, 101) = -10, (10, 51) = -1, (29, 1) = -7, (35, 101) = -54, (17, 1) = -65, (26, 51) = -39, (15, 101) = 32, (12, 151) = -1, (39, 101) = 87, (5, 1) = -10, (2, 51) = 29]);
convert(op(op(T)), set);
Matrix(40,151, %); 
 

One can prove that the first two equations for any nonnegative  u=z^(2^2011)   have the unique solution  x=y=u  (due to homogeneity, it is enough to check for  u=1) , which is clearly seen in the animation:

with(plots):
animate(implicitplot,[[x+y+x*y=u^2+2*u, x^4+y^4=2*u^4], x=-5..5,y=-5..5, color=[red,blue], axes=normal, gridrefine=5], u=0..3, frames=60, scaling=constrained);

Hence it is already easy to obtain the final solution   x=1, y=1, z=1

z1:=a+I*b: z2:=c+I*d:
evalc(z1*conjugate(z2)+conjugate(z1)*z2);
evalc(2*Re(z1*conjugate(z2)));

                                                       2*a*c+2*b*d
                                                       2*a*c+2*b*d

Two variants of the plotting: in a for loop for separate plots and all together on one plot:

for i from 30 by -5 to -30 do
plot([X(alpha1,i*Pi/180), Y(alpha1,i*Pi/180), alpha1=-30*Pi/180..30*Pi/180], view=[-1..1,-1..1]);
end do;   # The first option

plot([seq([X(alpha1,i*Pi/180), Y(alpha1,i*Pi/180), alpha1=-30*Pi/180..30*Pi/180], i=30..-30, -5)], view=[-1..1,-1..1]);   # The second option

You must save your document and then leave a link here using the bold green arrow in editor. Henceforth, do so as not to create unnecessary problems for those who respond to you.

I managed to recover your document using OCR system. The initial non-zero approximation was found in DirectSearch package, and then  fsolve  command was used to obtain the solution with high accuracy.

interface(rtablesize=infinity):
DirectSearch:-SolveEquations([E1,E2,E3], [a1=-5..5,a2=-5..5,a3=-5..5], AllSolutions);  # The best non-zero solution by DirectSearch

Digits:=20:
S1:=fsolve({E1,E2,E3}, {a1=3.53,a2=0.0075,a3=-1.312});  # The final result
eval([E1,E2,E3], S1);   # Check

For fsolve, I changed the sign of the initial  a1  and  a2, because the equations do not change when the signs change at  a1  and  a2 . 

System_solution.mw

First, each element of  S  should be converted to a list:

Matrix(convert~(S, list));

1. Two options:

f:=proc(x) if x<=1 then 1 else 2 fi;
end proc:

add(f(i), i=1..10);  # Or
sum('f'(i), i=1..10);

2. I did not understand at all what you are trying to do in the second example. For example, what does vars mean?

Or

Re(expr) assuming real;

Since your function at different intervals is given by different formulas, then to specify such a function, use  piecewise  command. See help on this command.

Addition - the simple example:

y[0]:=x^2;
y[1]:=1;
y[2]:=3-x;
y:=piecewise(x>=0 and x<1,y[0], x>=1 and x<2,y[1], x>=2 and x<=3,y[2]);
plot(y, x=0..3, scaling=constrained); 

You missed the arguments to the function  T_4  in the last line. Probably should be:

R11 := x -> (x+2*a_1)*(product(x+b[j], j = 1 .. k-n))/(x+a_1);
T_4 :=  (z_1, z_2, z_3, z_4) -> R11(z_1)*z_2*z_3*z_4;
%(z_1, z_2, z_3, z_4)*(z_1+a_1);
 

A simple recursive procedure called  FP  (fixed points)  finds the number of permutations of  n  elements that have exactly  m  fixed points. The procedure is written based on the formulas from the wiki (the link above):

FP:=proc(n,m)
option remember;
if m<0 or m>n then return 0 else
if n=0 and m=0 then return 1 else 
FP(n-1,m-1)+FP(n-1,m)*(n-1-m)+FP(n-1,m+1)*(m+1) 
fi; fi;
end proc:

Example of use:

seq(FP(10,i), i=0..10);
               1334961, 1334960, 667485, 222480, 55650, 11088, 1890, 240, 45, 0, 1

To achieve the same as on your sample graph, we can use some commands of  plottools  and  plots  packages:

restart:
with(plottools): with(plots):
h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
K1:=((4/h(z)^4)-(sin(alpha)/F)-h(z)^2+Nb*h(z)^4):
lambda:=Int(K1,z=0..1):
F:=0.3:
L1:=0.2:
d1:=0.2:
alpha:=Pi/6:
A:=plot( [seq(eval(lambda, Nb=j), j in [0.1,0.2,0.3])], delta2=0.02..0.1, axes=box, view=[0.02..0.1,1.6..2.66]):
B:=line([0.02,2.665],[0.10,2.665], thickness=0), seq(line([i,2.65],[i,2.665]), i=0.02..0.10,0.005):
C:=line([0.10,1.6],[0.10,2.665], thickness=0), seq(line([0.099,i],[0.1,i]), i=1.6..2.65,0.1):
display(A,B,C, axes=normal);

First, click on the capital letter  T , and then use Text  and  Math  buttons and the corresponding palettes if necessary.
Your example, done in this way:

For the numerical solution, we must specify the parameter values. I took  N=5, R=2 . Also we solve only  eq2  and  eq3  equations, since  eq1  is not a differential equation for  f(eta):

restart;
Digits := 15:
eq1:= diff(h(eta),eta)+2*f(eta):
eq2 := (1/4)*(diff(h(eta),eta))^2-g(eta)*g(eta)-(1/2)*h(eta)*diff(h(eta),eta,eta)+(1/2)*diff(h(eta),eta,eta,eta):
eq3 := -g(eta)*diff(h(eta),eta)+h(eta)*diff(g(eta),eta)-diff(g(eta),eta,eta): 
bc:=h(0)=0, D(h)(N)=0, D(h)(0)=-2, g(0)=R, g(N)=0:
N:=5; R:=2;
sol:=dsolve({eq2,eq3, bc}, numeric);
plots:-odeplot(sol,[[eta,-1/2*diff(h(eta),eta)], [eta,h(eta)], [eta,g(eta)]], eta=0..N, color=[red,blue,green], size=[600,400], scaling=constrained);

I suggest you another way of creating your animation through  plots:-animate  command. It's simpler, faster and allows you to easily adjust the animation speed programmatically simply by increasing the number of frames. In this method, you first write a procedure that creates one frame of the animation, and then just apply plots:-animate command.

restart;
OneFrame:=proc(R,r)
local sc, x_bar, bc, Pt, Rat, t1, t2;
uses plottools, plots;
Digits:=4;
sc := disk([0,2*R-r], r, color=white);
x_bar:=eval((R^2+R*r-r^2)/(R+r));
bc := disk([0,R], R, color=red);
Pt:=point([0,x_bar],color=black, symbol=solidcircle, symbolsize=10):
Rat:=R/r;
t1:=textplot([4*R/3,4*R/3,convert(x_bar, float)]);   
t2:=textplot([2*R,4*R/3,convert(Rat, float)]);
display(Pt,sc,t1,t2,bc, scaling=constrained, axes=none);
end proc:

plots:-animate(plots:-display, ['OneFrame'(20,r)], r=1..20, frames=120, paraminfo=false);  # Your animation of 120 frames