implicit real*8(a-h,o-z)
      common /lambda/al,eta,n
      dimension alv(20),d(20),xs(20)
c      open(1,file='feigenbaum.out',status='unknown')
      alv (1)=0.25
      alv(2)=0.75
      d(1)=0.5
      d(2)=0.25
      eta=0.1*2.5029
      dif=0.01
      do 10 n=1,13
      eta=eta/2.5029
      al=alv(n+1)+dif
      do 2 j=1,100
      x=ffp(dum)
      der=dn(x)
      t=1.+der
      if (t.le.O) goto 11
      al=al+dif
      tl=t
      all=al-dif
  2   continue
  11  tu=t
      alu=al
  12  al=.5*(all+alu)
      x=ffp(dum)
      der=dn(x)
      t=1.+der
      if (t*tl.gt.O) then
      all=al
      tl=t
      else
	alu=al
      tu=t
      endif
      if (alu-all.gt.10.**(-14)) goto 12
      alv(n+2)=al
      d(n+2)=abs(x-.5)
      alp=(d(n+1)-d(n))/(d(n+2)-d(n+1))
      del=(alv(n+1)-alv(n))/(alv(n+2)-alv(n+1))
      dif=.1*(alv(n+1)-alv(n))/del
c      write (1,22) n,del,alp,al,x
      print 22, n,del,alp,al,x
  22  format(i3,2f15.12,2f21.18)
  10  continue
      end


 
      function fn(y)
      implicit real*8(a-h,o-z)
      common /lambda/al,eta,n
      m=2**n
      bl=4.*al
      x=y
      do i=1,m
      x=bl*x*(1.-x)
      enddo
      fn=x
      return
      end

      function dn(y)
      implicit real*8(a-h,o-z)
      common /lambda/al,eta,n
      m=2**n
      bl=4.*al
      x=y
      d=1.
      do i=1,m
      d=d*bl*(1.-2.*x)
      x=bl*x*(1.-x)
      enddo
      dn=d
      return
      end

	function ffp(dum)
      implicit real*8(a-h,o-z)
      common /lambda/al,eta,n
      xl=.5-eta
      xu=.5+eta
      dd=fn(xl)-xl
      if (dd/(fn(xu)-xu).gt.O) then
      print *,n,'no solution'
      return
      endif
      do 1 i=1,40
      x=.5*(xl+xu)
      y=fn(x)
      if ((y-x)/dd.lt.0.) then
      xu=x	   
      else
      xl=x
      endif
  1   continue
      ffp=.5*(xl+xu)
      return
      end