!###########################################################
! page.75 表4.2(例題2)
! ニュートン・ラフソン法のプログラム
!
! f(x)=x**3+6*x**2+21*x+32
!###########################################################
PROGRAM prog14
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: KMAX=99
!-------------
x=0.d0
delta=1.d-10
!-------------
fx=func(x)
k=0
WRITE(*,"('x(',i2,')=',1PE15.8,', f(',i2,')=',1PE15.8)") k,x,k,fx
WRITE(*,*)


DO k=1,KMAX

xk=fnewton(x)
fx=func(xk)

WRITE(*,"('x(',i2,')=',1PE15.8,', f(',i2,')=',1PE15.8)") k,xk,k,fx

    IF(dabs(fx)<delta)GOTO 100

x=xk
    
END DO
WRITE(*,*) '解が求まりませんでした'
STOP
100 CONTINUE
WRITE(*,"('解は x=',1PE15.8)") xk
STOP
!###########################################################
CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(xx)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: func
   func=xx**3+6.d0*xx**2+21.d0*xx+32.d0
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION fnewton(xx)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: fnewton
fnewton=xx-(xx**3+6.d0*xx**2+21.d0*xx+32.d0)/(3.d0*xx**2+12.d0*xx+21.d0)
RETURN
END FUNCTION fnewton
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog14