c_fdjac_nn_sub Subroutine

subroutine c_fdjac_nn_sub(funcv,x,fjac,ml,mu,epsfcn)
  implicit none
  interface 
     subroutine funcv(n,x,y)
       integer                         :: n
       real(8),dimension(n),intent(in) :: x
       complex(8),dimension(n)         :: y
     end subroutine funcv
  end interface
  integer          ::  n
  real(8)          ::  x(:)
  complex(8)       ::  fvec(size(x))
  complex(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
  complex(8)       ::  wa1(size(x))
  complex(8)       ::  wa2(size(x))
  integer          :: i,j,k
  n=size(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
  call funcv(n,x,fvec)
  !  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
        call funcv(n,x,wa1)
        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
        call funcv(n,x,wa1)
        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)=dcmplx(0.d0,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 c_fdjac_nn_sub