*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     mufun.f
*
*     Copyright (C) 1991.
*
*     Marie Davidian and A. Ronald Gallant 
*     P.O. Box 5513, Raleigh NC 27650-5513, USA   
*
*     Permission to use, copy, modify, and distribute this
*     software and its documentation for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
      subroutine mufun(mode,ncnln,lvec,nrowj,needc,vec,c,cjac,nstate)   
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

      real*8 vec,c,cjac
      real*8 work
      real*8 R0
      real*8 Eyl,dEyl

      integer*4 mode,ncnln,lvec,nrowj,needc,nstate
      integer*4 lwork,mwork
      integer*4 lambda
      integer*4 i,j,ii,jj

      parameter (mwork = 500)

      dimension needc(*)
      dimension vec(lvec),c(*),cjac(nrowj,*)
      dimension R0((mM*(mM+1))/2)
      dimension lambda(mM)
      dimension dEyl(mtheta)
      dimension work(mwork)


      if (ncnln.eq.0) return

      lwork=2*ntermz
      if (lwork.gt.mwork) then
        write (3,3906) lwork
        stop
      endif

*     Copy parameters NPSOL gets to see to vsav.

      j=0
      do i=1,lvsav
        if (ivec(i).eq.0) then
          j=j+1
          vsav(i)=vec(j)
        end if
      end do 

*     Make tau and theta agree with vsav.

      do i=1,ltau
        tau(i)=vsav(i)
      end do

      do i=1,ltheta
        theta(i)=vsav(ltau+i)
      end do

*     This is an identity to pass to EYLX instead of R.

      do j=1,M
      do i=1,j
        if (i.eq.j) then
          R0(j*(j-1)/2+i)=1.d0
        else
          R0(j*(j-1)/2+i)=0.d0
        end if
      end do
      end do

      do i=1,M

        if (needc(i).gt.0) then

          do jj=1,M
            lambda(jj)=0
          end do
          lambda(i)=1

          call eylx(lambda,z,x,Eyl,dEyl,
     &      theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &      M,L,ntermz,ntermx,Kz,Kx,work)

          c(i) = Eyl

          jj=0
          do ii=1,lvsav
            if (ivec(ii).eq.0) then
              jj=jj+1
              if ((ltau+ibA+1.le.ii).and.(ii.le.ltau+ibA+lA)) then
                 cjac(i,jj)=dEyl(ibA+ii-ltau)
              else
                cjac(i,jj)=0.d0
              end if
            end if
          end do 

        end if
      end do

      return
3906  format(1x,'error, muzero, increase mwork to at least',i5)
      end
