syms a X f
f=a*X-X^3                                     % given function, x_{n+1}=f(x_n)
fixed_pts=solve(f-X,X)                        % solve f(X)-X=0
deriv=diff(f,X)                               % differentiate f
deriv_fixed_pts=subs(deriv,X,fixed_pts)       % evaluate f at the fixed points
two_cycle=collect(a*f-f^3-X,X)                % set up f(f(X))-X
pretty(two_cycle)
factors=factor(two_cycle)                     % factorize f(f(X))-X
pretty(factors)
factors=simplify(factors/(X^3+(1-a)*X))       % remove known simple fixed points
roots2=solve(factors,X)                       % find the six possible two-cycles
solve(two_cycle,X);                           % alternative approach: also gives
                                              % fixed points
syms X1 X2 f1 f2
f1=a*X1-X1^3;
f2=a*X2-X2^3;
stable=diff(f1,X1)*diff(f2,X2)                % stability criterion for 2-cycle
crit1=simplify(subs(stable,[X1;X2],[roots2(1);roots2(4)])) % evaluate stability
                                                           % criterion for 1st pair
crit2=simplify(subs(stable,[X1;X2],[roots2(2);roots2(3)])) % evaluate stability
                                                           % criterion for 2nd pair
crit3=simplify(subs(stable,[X1;X2],[roots2(5);roots2(6)])) % evaluate stability
                                                           % criterion for 3rd pair