!###########################################################
! page.87-88 表4.5(例題6)
! マラー法のプログラム
!
! f(x)=x**3+6.d0*x**2+21.d0*x+32.d0
!###########################################################
PROGRAM prog18
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: KMAX=100

!---------------
x0=-3.d0
x1= 0.d0
x2=-2.d0

delta=1.d-8
!---------------
y0=func(x0)
y1=func(x1)
y2=func(x2)

h=(x2-x1)
r=(x2-x1)/(x1-x0)

DO k=3,KMAX

g=(1.d0+2.d0*r)*(y2-y1)-r**2*(y1-y0)
D=g**2-4.d0*y2*(1.d0+r)*r*((y2-y1)-r*(y1-y0))
  
  IF(dabs(g+dsqrt(D))>=dabs(g-dsqrt(D)))THEN
  gd=g+dsqrt(D)
  ELSE
  gd=g-dsqrt(D)
  END IF

r=-2.d0*y2*(1.d0+r)/gD

!------------------
x3=x2+h*r
y3=func(x3)
h=h*r
!------------------
WRITE(*,"('k=',i2,', x=',1PE16.8,', f(x)=',1PE16.8)") k,x3,y3
IF(dabs(y3)<delta) GOTO 100

x0=x1 ; y0=y1
x1=x2 ; y1=y2
x2=x3 ; y2=y3

END DO
WRITE(*,*) '解が求まりませんでした'
STOP

100 CONTINUE
WRITE(*,*)
WRITE(*,*) "マラー法で求めた実数解は"
WRITE(*,"('x=',1PE15.8)") x3

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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog18
!###########################################################