fdjac2 Subroutine

    subroutine fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, Iflag, Epsfcn, Wa)

        implicit none

        procedure(func2) :: fcn !! the user-supplied subroutine which
                            !! calculates the functions.
        integer, intent(in) :: m !! a positive integer input variable set to the number
                            !! of functions.
        integer, intent(in) :: n !! a positive integer input variable set to the number
                            !! of variables. n must not exceed m.
        integer, intent(in) :: Ldfjac !! a positive integer input variable not less than m
                                 !! which specifies the leading dimension of the array fjac.
        integer, intent(inout) :: Iflag !! an integer variable which can be used to terminate
                                 !! the execution of fdjac2. see description of [[func2]].
        real(wp), intent(in) :: Epsfcn !! an input variable used in determining a suitable
                                  !! step length for the forward-difference approximation. this
                                  !! approximation assumes that the relative errors in the
                                  !! functions are of the order of epsfcn. if epsfcn is less
                                  !! than the machine precision, it is assumed that the relative
                                  !! errors in the functions are of the order of the machine
                                  !! precision.
        real(wp), intent(inout) :: x(n) !! an input array of length n.
        real(wp), intent(in) :: Fvec(m) !! an input array of length m which must contain the
                                   !! functions evaluated at x.
        real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output m by n array which contains the
                                            !! approximation to the jacobian matrix evaluated at x.
        real(wp), intent(inout) :: Wa(m) !! a work array of length m.

        integer :: i, j
        real(wp) :: eps, h, temp

        eps = sqrt(max(Epsfcn, epsmch))
        do j = 1, n
            temp = x(j)
            h = eps*abs(temp)
            if (h == zero) h = eps
            x(j) = temp + h
            call fcn(m, n, x, Wa, Iflag)
            if (Iflag < 0) return
            x(j) = temp
            do i = 1, m
                Fjac(i, j) = (Wa(i) - Fvec(i))/h
            end do
        end do

    end subroutine fdjac2