subroutine fdjac_nn_func(funcv,x,fjac,ml,mu,epsfcn)
interface
function funcv(x)
real(8),dimension(:),intent(in) :: x
real(8),dimension(size(x)) :: funcv
end function funcv
end interface
integer :: n
real(8),intent(in) :: x(:)
real(8) :: x_(size(x))
real(8) :: fvec(size(x))
real(8) :: fjac(size(x),size(x))
integer,optional :: ml, mu
real(8),optional :: epsfcn
integer :: ml_,mu_,msum
real(8) :: eps,eps_
real(8) :: epsmch
real(8) :: h,temp
real(8) :: wa1(size(x))
real(8) :: wa2(size(x))
integer :: i,j,k
n=size(x)
x_ = x
ml_ = n-1 ; if(present(ml))ml_=ml
mu_ = n-1 ; if(present(mu))mu_=mu
eps_= 0.d0; if(present(epsfcn))eps_=epsfcn
epsmch = epsilon(epsmch)
eps = sqrt(max(eps_,epsmch))
msum = ml_ + mu_ + 1
! Evaluate the function
fvec = funcv(x_)
! Computation of dense approximate jacobian.
if(n <= msum)then
do j=1,n
temp = x_(j)
h = eps*abs(temp)
if(h==0.d0) h = eps
x_(j) = temp + h
wa1 = funcv(x_)
x_(j) = temp
fjac(1:n,j) = ( wa1(1:n) - fvec(1:n) )/h
enddo
else
! 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==0.d0)h = eps
x_(j) = wa2(j) + h
end do
wa1 = funcv(x_)
do j=k,n,msum
x_(j) = wa2(j)
h = eps*abs(wa2(j))
if(h==0.d0)h = eps
enddo
fjac(1:n,j)=0.d0
do i=1,n
if( (j-mu_<=i).AND.(i<=j+ml_) )then
fjac(i,j) = ( wa1(i) - fvec(i) )/h
end if
end do
end do
end if
return
end subroutine fdjac_nn_func