> restart;
# Period three saddle-node bifurcation for the Logistic Map
> f := x -> r*x*(1-x);
# We note that f^3(x) crosses the diagonal somewhere between r=3.82 and
# 3.84.
> plot([x,subs(r=3.82,f(f(f(x)))),subs(r=3.84,f(f(f(x))))],x=0..1);
# To find the exact value, we look at the equation for the fixed points
> p3a := factor(x -f(f(f(x))));

# Note: this has as roots x=0 and x=(1-r)/r, the fixed points. So we
# divide them out
> p3b := p3a/(x*(r*x+1-r));
# Now we find the equation for the stability multiplier
# First differentiate f
> fp := unapply(diff(f(x),x),x);
# Now the multiplier is the product of the derivatives along the orbit

> x1 := f(x);
> x2 := f(x1);
> lambda := collect(fp(x)*fp(x1)*fp(x2),r); 
# Here is the big trick. The "Resultant" of two polynomials with
# respect to a variable x gives a polynomial that eliminates x, and is
# zero only when both polynomials vanish.
# Here we eliminate the point x in favor of the parameter r. Note that
# we what to find when lambda = 1, so we use lambda-1
#  
> eq:= resultant(lambda-1,p3b,x);
> factor(eq);
# The resultant, amazingly, has a number of simple factors. One of
# these must be zero in order that lambda = 1. We can solve each of them
> solve(r^2-2*r-7,r);
> solve(r^2-5*r+7,r);
> solve(r^2+r+1,r);
> evalf(1+sqrt(8));
# Note that only the first one has real roots. Thus we get a
# saddle-node bifurcation at r = 1 ± sqrt(8). To find the value of x
# that this happens we solve p3b:
> Digits := 10;
> fsolve(subs(r=1+sqrt(8),p3b),x=0.15,fulldigits);
> evalf(subs(r=1+sqrt(8),x=.15992881838665496465,p3b));
> evalf(subs(r=1+sqrt(8),x=.15992881850585781331,p3b));
# Curiously, Maple thinks there are two nearby roots.
> plot(subs(r=1+sqrt(8),p3b),x=0.1599..0.160,y=-1.0e-7..1.0e-7);
# It appears that this function has a unique zero, but it is hard to
# tell. Changing the number of "Digits" above, changes the two answers.
# So Maple is making some numrical error.

# Similarly, we could solve for the period doubling point of the period
# 3 orbit


> eqpd:= resultant(lambda+1,p3b,x);
> factor(eqpd);
# This doesn't have any simple analytical solutions. We can get a
# numerical value.
> fsolve(eqpd/r^42, r=3.8..3.95);