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

@brian bovril

I do not know the closed formula for the number of partitions of a set into subsets of different lengths, but if the sizes of subsets are known, then the number of ways to split a set into subsets of specified sizes can be found using the following procedure:

NumbPart:=proc(Q::list(posint))  # Procedure finds the number of all partitions of a set into subsets of the specific size given  Q
local L, T, P, n, S, N;
uses ListTools;
L:=convert(Q, set);
T:=[seq([L[i], Occurrences(L[i], Q)], i=1..nops(L))];
P:=convert(T, set);
n:=nops(P);  N:=add(P[i,2], i=1..n);
S:=add(P[i,1]*P[i,2], i=1..n)!/mul(P[i,1]!^P[i,2], i=1..n)/mul(P[i,2]!, i=1..n);
end proc:

Examples of use:

NumbPart([3,3,3]); 
NumbPart([2,3,4]);
NumbPart([10,10,10,10,10]);
                                                 280
                                                1260
                        402789797982510165934296910320

If it is necessary to count the total number of partitions of a certain set (it's size is  n) into subsets of different length, then first, using  Partition  procedure (see above in this thread), we find possible variants of partitioning the number n, then select of them the partitions of different lengths, and then applying  NumbPart procedure, we obtain the final result.

Example: Find the number of partitions of a set of 100 elements into 3 subsets of different lengths.

Partition(100, 3):
select(s->nops(s)=nops({s[]}), %);
`+`(NumbPart~(%)[]);

                                    73329120859641139777709300892745606087338612215

g:=unapply(int(exp(z^2), z = 0 .. x), x);   # The initial function
h:=unapply(solve(y=g(x), x, explicit), y);   # The inverse function

g(2);   # Examples of use
h(%);
evalf(%);

Unfortunately this does not work for more complex cases. Here is a workaround (in Maple 2017.3):

restart;
`union`(solve({x^2-3*x+2>0, x>=0, x<10}));     
`union`(solve(x^2-3*x+2>0 and x>=0 and x<10));

It seems that the use of this syntax ( and  and/or  or  in  solve  command) is not documented in Maple, but it is useful in some cases.