      PROGRAM bifurcate
*     plot values of x for different values of r
      IMPLICIT NONE
      REAL*8 dr, r, rmax, x
      REAL r0
      INTEGER IR_, ir, nplot, ntransient, nvalues
      CALL GWopen(IR_, 0)
      CALL initial(x,r,rmax,nvalues,dr,ntransient,nplot)
      r0 = REAL(r)
      DO ir = 0, nvalues
         CALL output(x,r,ntransient,nplot,r0)
         r = r + dr
      END DO
*     maximum value of r done separately to avoid r > 1
      CALL output(x,rmax,ntransient,nplot,r0)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(x0,r0,rmax,nvalues,dr,ntransient,nplot)
      IMPLICIT NONE
      REAL*8 x0, r0, rmax, dr
      REAL mx, xmax
      INTEGER IR_, nvalues, ntransient, nplot
      WRITE(*,'(A,$)') 'initial value of control parameter r = '
      READ(*,*) r0
*     important that r not be greater than 1
      WRITE(*,'(A,$)') 'maximum value of r = '
      READ(*,*) rmax
*     suggest dr <= 0.01
      WRITE(*,'(A,$)') 'incremental change of r = '
      READ(*,*) dr
      WRITE(*,'(A,$)') 'number of iterations not plotted = '
      READ(*,*) ntransient
      WRITE(*,'(A,$)') 'number of iterations plotted = '
      READ(*,*) nplot
      nvalues = int((rmax - r0)/dr)         ! number of r values plotted
      x0 = 0.5D0                            ! initial value
      xmax = 1.0                            ! maximum value of x
      mx = 0.05*xmax                        ! margin
      CALL GWsetogn(IR_, 0, 1)
      CALL GWindow(IR_, REAL(r0-dr), -mx, REAL(rmax+dr), xmax + mx)
      CALL GWrect(IR_, REAL(r0), 0.0, REAL(rmax), 1.0)
      END
C
      SUBROUTINE output(x,r,ntransient,nplot,r0)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 x, r, f
      REAL r0
      INTEGER IR_, i, ntransient, nplot
      CALL GWsetogn(IR_, 0, -2)
      WRITE(Stmp1_,*) 'r = ', r
      CALL GWputtxt(IR_, r0, 1.0, Stmp1_)
      CALL GWflush(IR_, -2)
      CALL GWsetogn(IR_, 0, 1)
      DO i = 1, ntransient
         x = f(x,r)
      END DO
      DO i = 1, 0.5*nplot
         x = f(x,r)
*        show different x-values for given value of r
         CALL GWsetpxl(IR_, REAL(r), REAL(x), 13)
      END DO
*     change color to see if values of x have converged
      DO i = (0.5*nplot + 1), nplot
         x = f(x,r)
         CALL GWsetpxl(IR_, REAL(r), REAL(x), 16)
      END DO
      END
C
      FUNCTION f(x,r)
      IMPLICIT NONE
      REAL*8 x, r, f
      f = 4*r*x*(1 - x)
      END
