secant_solve Subroutine

  subroutine secant_solve(f,x1,x2,solution,error,residual,opt_tol)
      ! Solve an equation f(x)=0
      ! Only requires one function evaluation per iteration (plus two to start with)
      ! Does not require any derivatives.
      ! https://en.wikipedia.org/wiki/Secant_method
      ! Requires two initial values, x0 and x1, which should ideally be chosen to lie close to the root.
      !external:: f

      interface
        function f(x)
#ifndef dp
          use, intrinsic :: iso_fortran_env, only: dp=>real64
#endif
          real(dp), intent(in) :: x
          real(dp) :: f
        end function f
      end interface

      external :: f
      real(dp), intent(out) ::solution, residual
      real(dp), intent(in) ::x1,x2
      real(dp), intent(in), optional ::opt_tol
      real(dp),dimension(20) ::x
      integer :: i
      logical, intent(out)::error
      real(dp)::mean, tol,fximinus1, fximinus2
      error = .FALSE.
      tol=0.001d0; if (present(opt_tol)) tol=opt_tol

      x(1)=x1
      x(2)=x2
      mean = (x1+x2)/2
      ! Calculate the first two values before the loop
      fximinus1 = f(x(2))
      fximinus2 = f(x(1))
      !write(*,*)"x(1)",x(1),"x(2)",x(2),"fximinus1",fximinus1,"fximinus2",fximinus2
      do i=3,20
          x(i) = x(i-1) - fximinus1 * ((x(i-1) - x(i-2)) / (fximinus1 - fximinus2))
          ! Store values for the *next iteration
          fximinus2 = fximinus1
          fximinus1 = f(x(i))
          residual = fximinus1
          !write(*,*)"x(i)", x(i), "  residual",residual,"  fximinus1",fximinus1, "  fximinus2",fximinus2
          ! test for convergence
          if(abs(residual) < tol) then
              solution = x(i)
              return
          endif
      end do
      ! Convergence not achieved.  Return the best value and a flag.
      error = .TRUE.
      solution = x(i-1)
      write(*,*)"Secant solver not converged.  solution", solution, "  residual",residual
      !stop 1

  end subroutine secant_solve