PROGRAM period
! find fixed point of f(x) iterated p times
DECLARE DEF f
CALL initial(r,p,epsilon,xleft,xright,gleft,gright)
DO
   CALL bisection(r,p,xleft,xright,gleft,gright)
LOOP until abs(xleft - xright) < epsilon
LET x = 0.5*(xleft + xright)
PRINT "explicit demonstration of period"; p; "behavior"
PRINT
PRINT 0,x                    ! result
FOR i = 1 to 2*p + 1
    LET x = f(x,r,1)
    PRINT i,x
NEXT i
END

SUB initial(r,p,epsilon,xleft,xright,gleft,gright)
    DECLARE DEF f
    INPUT prompt "control parameter r = ": r
    INPUT prompt "period = " : p
    INPUT prompt "desired precision = ": epsilon
    LET done$ = "no"
    DO while done$ = "no"    ! do until zero between xleft and xright
       INPUT prompt "guess for xleft = ": xleft
       INPUT prompt "guess for xright = ": xright
       LET gleft = f(xleft,r,p) - xleft
       LET gright = f(xright,r,p) - xright
       IF gleft*gright < 0  then
          LET done$ = "yes"
       ELSE
          PRINT "range does not necessarily enclose a root"
       END IF
    LOOP
END SUB

SUB bisection(r,p,xleft,xright,gleft,gright)
    DECLARE DEF f
    ! midpoint between xleft and xright
    LET xmid = 0.5*(xleft + xright)
    LET gmid = f(xmid,r,p) - xmid
    IF gmid*gleft > 0 then
       LET xleft = xmid      ! change xleft
       LET gleft = gmid
    ELSE
       LET xright = xmid     ! change xright
       LET gright = gmid
    END IF
END SUB

DEF f(x,r,p)                 ! f defined by recursive procedure
    IF p > 1 then
       LET y = f(x,r,p-1)
       LET f = 4*r*y*(1-y)
    ELSE
       LET f = 4*r*x*(1-x)
    END IF
END DEF