!===========================================================
!Newton Raphson of function of one variable
!===========================================================
program NewtonRaphson
 implicit none
 real(8)            :: r
 real(8), parameter :: eps  = 1D-12
 integer, parameter :: NMAX = 100
 real(8) :: x0, x1, err, g, gp
 integer :: i
 integer :: k
 real(8) :: f, dfdx

 print *, 'Enter k,r, x0: '
 read  *, k, r, x0

 err = 1.0D0
 print *,'iter           x                        error    '
 print *,'-------------------------------------------------'
 print *, 0,x0,err
 do i=1,NMAX
!value of function
  g   = x0    - f(k,r,x0)
!value of the derivative
  gp  = 1.0D0 - dfdx(k,r,x0)
  x1  = x0 - g/gp
  err = ABS(x1-x0)
  print *,i,x1,err
  if(err .lt. eps) exit
  x0 = x1
 enddo
 print *,'-------------------------------------------------'
 print *,'k=',k,' r= ',r,' x* = ',x1,' |fp(x*)|=',abs(dfdx(k,r,x1))
end program NewtonRaphson
!-----------------------------------------------------
real(8) function f(k,r,x)
  implicit none
  integer :: k,i
  real(8) :: r,x,x0

  x0=x
  do i=1,k
     x0 = r*x0*(1.0D0-x0)
  enddo
  f = x0
end function f
!-----------------------------------------------------
real(8) function dfdx(k,r,x)
  implicit none
  integer :: k
  real(8) :: x,r,eps
  real(8) :: f
  
  eps  = 0.001D0 *x
  dfdx = (f(k,r,x+eps)-f(k,r,x-eps))/(2.0D0*eps)

end function dfdx
