!===========================================================
!Newton Raphson for bifurcation points: File nr3.f90 (in NotesHandwritten/Lecture07/)
!===========================================================
program NewtonRaphson
 implicit none
 integer , parameter :: dp   = 8
 real(dp), parameter :: eps  = 1.0e-12_dp
 integer , parameter :: NMAX = 100
 real(dp)            :: err  = 1.0_dp
 real(dp)            :: r
 real(dp)            :: x0, x1, g, gp
 integer             :: i
 integer             :: k
 real(dp), parameter :: ONE=1.0_dp,TWO=2.0_dp

 print *, 'Enter k, r, x0: '
 read  *,        k, r, x0
 print *,'iter           x                        error    '
 print *,'-------------------------------------------------'
 print *, 0,x0,err
 do i=1,NMAX
!value of function         g (x) = x - r*x*(1-x)
  g   = x0    - f   (x0)
!value of the derivative   g'(x) = 1 - r + 2*r*x
  gp  = ONE   - dfdx(x0)
  x1  = x0    - g/gp
  err = ABS(x1-x0)
  print *,i,x1,err
  if(err < eps) exit
  x0  = x1
 enddo
 print *,'-------------------------------------------------'
 print *,'k=',k,' r= ',r,' x* = ',x1,' |fp(x*)|=',abs(dfdx(x1))
!-----------------------------------------------------
contains
!-----------------------------------------------------
 function f(x)
  real(dp) :: f
  real(dp) :: x,x0
  integer  :: i
  x0   = x
  do i = 1,k
   x0  = r*x0*(ONE-x0)
  enddo
  f    = x0
 end function f
!-----------------------------------------------------
 function dfdx(x)
  real(dp) :: dfdx
  real(dp) :: x,eps
  
  eps  = 0.001_dp *x

  dfdx = (f(x+eps)-f(x-eps))/(TWO*eps)

 end function dfdx
!-----------------------------------------------------
end program NewtonRaphson
