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.

restart;
A:=(-cos(theta) + 1)*z^2  + cos(theta):
local cos: 
B:=op([1,2], A)*(1-cos(theta))+op(2, A):
B;

You are plotting a curve in space, not in a plane. Therefore, replace the last line of code with

plots:-spacecurve([x__1, x__2, x__3], t2 = 0 .. (1/3)*tau, color=red);

azido_displacement_new.mw

Replace the last line of your code with two lines:

V2:=x+I*y: V3:=u+I*v:
solve({eq1, eq2}, {x,y,u,v}) assuming real;
                                        {u = -1.*I*v, v = v, x = -1.*I*y, y = y}

The result obtained means that there are an infinite number of solutions depending on two parameters  y  and  v .

Use  plots:-contourplot  command instead. See help on this command.

PS. If you still want to use  listcontplot  command, then use a linear transformation  [1..m, 1..n]->[x1..x2, y1,,y2], which is given by the obvious formulas

F:=(x,y)->[x1+(x2-x1)/(m-1)*(x-1), y1+(y2-y1)/(n-1)*(y-1)]

and then  plottools:-transform command.

Edit.

For not too large numbers, I think  length  command is the best option. An alternative to length is ilog10  command:

length(12345678910);
ilog10(12345678910)+1;
                                                        11
                                                        11

Moreover, length command probably uses  ilog10  in its code. However, both commands fail for very large numbers:

length(12345^(123456789));
ilog10(12345^(123456789))+1;
        Error, numeric exception: overflow
        Error, numeric exception: overflow

The workaround is very simple - you just need to use the well-known property of logarithms:

floor(123456789*evalf[20](log10(12345))+1);
                                                           505122353

Y := {Theta[i, j, 2]+2*Omega[i, j, 2]+4*GAMMA[i, j, 2], Omega[i, j, 1]+3*GAMMA[i, j, 1]+Omega[i, j, 2]+3*GAMMA[i, j, 2], -Theta[i, j, 1]+2*GAMMA[i, j, 1]+Omega[i, j, 2]+3*GAMMA[i, j, 2]};
nops(Y);