fdjac1 Subroutine

    subroutine fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, Iflag, Ml, Mu, Epsfcn, Wa1, Wa2)

        implicit none

        procedure(func) :: fcn !! the user-supplied subroutine which
                            !! calculates the functions.
        integer, intent(in) :: n !! a positive integer input variable set to the number
                            !! of functions and variables.
        integer, intent(in) :: Ldfjac !! a positive integer input variable not less than n
                                 !! 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 fdjac1. see description of [[func]].
        integer, intent(in) :: Ml !! a nonnegative integer input variable which specifies
                             !! the number of subdiagonals within the band of the
                             !! jacobian matrix. if the jacobian is not banded, set
                             !! ml to at least n - 1.
        integer, intent(in) :: Mu !! a nonnegative integer input variable which specifies
                             !! the number of superdiagonals within the band of the
                             !! jacobian matrix. if the jacobian is not banded, set
                             !! mu to at least n - 1.
        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(n) !! an input array of length n which must contain the
                                   !! functions evaluated at x.
        real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output n by n array which contains the
                                            !! approximation to the jacobian matrix evaluated at x.
        real(wp), intent(inout) :: Wa1(n) !! work array of length n.
        real(wp), intent(inout) :: Wa2(n) !! work array of length n. if ml + mu + 1 is at
                                     !! least n, then the jacobian is considered dense, and wa2 is
                                     !! not referenced.

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

        eps = sqrt(max(Epsfcn, epsmch))
        msum = Ml + Mu + 1
        if (msum < n) then
            ! computation of banded approximate jacobian.
            do k = 1, msum
                do j = k, n, msum
                    Wa2(j) = x(j)
                    h = eps*abs(Wa2(j))
                    if (h == zero) h = eps
                    x(j) = Wa2(j) + h
                end do
                call fcn(n, x, Wa1, Iflag)
                if (Iflag < 0) return
                do j = k, n, msum
                    x(j) = Wa2(j)
                    h = eps*abs(Wa2(j))
                    if (h == zero) h = eps
                    do i = 1, n
                        Fjac(i, j) = zero
                        if (i >= j - Mu .and. i <= j + Ml) Fjac(i, j) = (Wa1(i) - Fvec(i))/h
                    end do
                end do
            end do
        else
            ! computation of dense approximate jacobian.
            do j = 1, n
                temp = x(j)
                h = eps*abs(temp)
                if (h == zero) h = eps
                x(j) = temp + h
                call fcn(n, x, Wa1, Iflag)
                if (Iflag < 0) return
                x(j) = temp
                do i = 1, n
                    Fjac(i, j) = (Wa1(i) - Fvec(i))/h
                end do
            end do
        end if

    end subroutine fdjac1