Program FindRoot
 Implicit None
 Real :: a, b, x, ft_f
 Integer :: maxit,id

 Do
  Print *, " Bisection(1), Newton(2), Stop(0) "
  Read *, id
  If (id==0) Then
     Stop
  ElseIf (id==1) Then
     Print *, " Input bound (a, b) = "
     Read *, a,b
  Else
     Print *, " Input an initial guess x0 = "
     Read *, a
  EndIf
     
!  Print *, " Type the # of iterations maxit = "
!  Read *, maxit
!  If (maxit==0) Stop
  maxit = 200

  If (id == 2) then
    Call Newton(x,a,maxit)
  Else
    Call Bisection(x,a,b,maxit)
  EndIf

  print *, "------------------------------------------"
  If (id==2) then
    print *, " The Result by Newton Iteration "
  Else
    print *, " The Result by Bisection Iteration "
  EndIf
  write (*,"(A20,pE23.16)") " The root is ", x
  write (*,"(A20,pE23.16)") " The value is ", ft_f(x)
  write (*,"(A20,I4)")     " The iteration is ", maxit 
  print *, "------------------------------------------"
 EndDo

End Program FindRoot
  
Subroutine Bisection(x,a,b,maxit)
  Implicit None
  Real, Parameter :: tol=0.1**15
  Real :: x,a,b,v
  Integer :: maxit, i
  Real :: ft_f

  If (ft_f(a)*ft_f(b)>0) then 
    Print *, " Must be f(a)*f(b) < 0 "
    Return
  EndIf
  If (ft_f(a)==0) then ; x = a ; Return ; EndIf
  If (ft_f(b)==0) then ; x = b ; Return ; EndIf

  If (ft_f(a)>0) then ; x=a; a=b; b=x; EndIf
  
  Do i=1,maxit
    x = (a+b)/2
    v = ft_f(x)
    Write(*,"(I4,A,pE23.16,A,pE23.16)") i,"-th x = ",x," value = ",ft_f(x) 
    If (v<-tol) then ; a=x
      ElseIf (v>tol) then ; b=x
      Else ; maxit = i ; Return
    EndIf
  EndDo

End Subroutine Bisection
  
Subroutine Newton(x,a,maxit)
  Implicit None
  Real, Parameter :: tol=0.1**15
  Real :: x,a,v
  Integer :: maxit, i
  Real :: ft_f, ft_fp ! f(x), f'(x)

    Write(*,"(I4,A,pE23.16,A,pE23.16)") 0,"-th x = ",a," value = ",ft_f(a) 
  Do i=1,maxit
    If (ft_fp(a)==0) Exit
    x = a - ft_f(a)/ft_fp(a)
    a = x
    Write(*,"(I4,A,pE23.16,A,pE23.16)") i,"-th x = ",x," value = ",ft_f(x) 
    If (abs(ft_f(x))<tol) Exit
  EndDo
    maxit = i

End Subroutine Newton
  
Real Function ft_f(x)
  Implicit None
  Real x

  ft_f = x**2+2*x+sin(x)*exp(-x)
End Function ft_f

Real Function ft_fp(x)
  Implicit None
  Real x

  ft_fp = 2*x+2+cos(x)*exp(-x)-sin(x)*exp(-x)
End Function ft_fp