See the corrected file Code1_new1.mws

The last line of code does not plot anything, because  g1(t)  is not defined.

Edit - the file was replaced.

In fact, your system is very simple: one equation is linear, and the second is quadratic with respect to the indicated unknowns. Therefore, all solutions are easily found by  solve  command (I took the same parameter values that in acer's solution):

restart;
local gamma:
eq1 := gamma=xx*sinh(xx)+xi*psi*cosh(xx):
eq2 := psi=-gamma^2+beta+sqrt(xx/w)*coth(sqrt(xx/w))-xx:
params:=[beta=1.0,w=2.0,xi=3.0,xx=3.0]:
solve(eval([eq1,eq2],params),{gamma,psi});

  {gamma = 0.6554134802, psi = -0.9733544660}, {gamma = -0.6885227893, psi = -1.017851267}

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: