PROGRAM Newtons_Method
!-----------------------------------------------------------------------  
! Program to find an approximate root of a function F using Newton's   
! method.  Variables used are:
!   F, FPrime     : the function and its derivative (internal functions)
!   OldApprox     : previous approximation (initially the first one)
!   FDeriv_Old    : value of the derivative of F at OldApprox
!   NewApprox     : the new approximation
!   F_Value       : value of F at an approximation
!   Epsilon       : repetition stops when ABS(F_Value) is less than
!                   Epsilon
!   NumIterations : limit on number of iterations
!   N             : number of iterations
!
! Input:  Epsilon,NumIterations, and OldApprox
! Output: Iteration number N, the Nth approximation, and the value of
!           F at that approximation, or an error message indicating 
!           that the method fails
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: NumIterations, N
  REAL :: OldApprox, FDeriv_Old, NewApprox, Epsilon, F_Value

  ! Get termination values Epsilon and NumIterations and 
  ! initial approximation

  PRINT *, "Enter epsilon , limit on # of iterations,"
  PRINT *, " and the initial approximation:"
  READ *, Epsilon, NumIterations, OldApprox
  PRINT *

  ! Initialize F_Value and N; print headings and initial values

  F_Value = F(OldApprox)
  N = 0
  PRINT *, "  N      X(N)       F(X(N))"
  PRINT *, "============================="
  PRINT 10, 0, OldApprox, F_Value
  10 FORMAT(1X, I3, F11.5, E14.5)

! Iterate using Newton's method while ABS(F_Value) is greater
! than or equal to Epsilon and N has not reached NumIterations

  DO
     IF ((ABS(F_Value) < Epsilon) .OR. (N > NumIterations)) EXIT
     ! If a termination condition met, stop generating approximations

     ! Otherwise continue with the following
     N = N + 1
     FDeriv_Old = FPrime(OldApprox)

     ! Terminate if the derivative is 0 at some approximation
     IF (FDeriv_Old == 0) THEN
        PRINT *, "Newton's method fails -- derivative = 0"
        EXIT
     END IF

     ! Generate a new approximation
     NewApprox = OldApprox - (F_Value / FDeriv_Old)
     F_Value = F(NewApprox)
     PRINT 10, N, NewApprox, F_Value
     OldApprox = NewApprox

  END DO

CONTAINS

  !-F(X)--------------------------------------
  ! Function for which a root is being found  
  !-------------------------------------------
  
  FUNCTION F(X)
  
    REAL :: F
    REAL, INTENT (IN) :: X 
  
    F = X**3 + X - 5.0
  
  END FUNCTION F
  
  
  !-FPrime(X)---------------------------------
  ! The derivative of the function F
  !-------------------------------------------
  
  FUNCTION FPrime(X)
  
    REAL :: FPrime
    REAL, INTENT (IN) :: X 
  
    FPRIME = 3.0*X**2 + 1.0
  
  END FUNCTION FPrime
  
END PROGRAM Newtons_Method