restart:
with(plots):H := x -> r*x/(1+x)^b;
dH :=x -> r*(1+x-b*x)/(1+x)^(1+b);r := 10.0;
b := 5;
n := 16:
x0 := 1/(b-1);
xmax := 1.2*H(x0);
p := evalf(r^(1/b)-1);
evalf(dH(p));
xc := evalf(1/(b-1));
xc1 := H(xc);
xc2 := H(xc1);pic := {plot({H(x),x}, x=0..xmax,color = black,thickness=2)}:
for i from 1 to n do
x1:= H(x0):
pic := pic union {plot([[x0,x0],[x0,x1],[x1,x1]], color = black,thickness=2)}:
x0:=x1:
od:
display(pic);a := evalf(3^(3/2));g := z -> -a*z^18 + 4+a^3 + 12*a^2*z + 12*a*z^2 +4*z ;g(a);
plot(g(z), z= a..2*a);
r := 40:
b := 5:
n := 50:
x0 := evalf(1/(b-1));
for i from 1 to n do
x0:= H(x0);
od;