c  ---------------     please, cut here    --------------------------


c
c  example for the call of subroutine limex
c  integrator for differential-algebraic systems
c
      implicit real*8(a-h,o-z)
      double precision y(3),rw(150)
      integer iw(30),ijob(13)
      external fcn,lsc,jac,jdum
      n=3
      nzc=1
      nzv=1
      t=0.d0
      y(1)=-1.d0
      y(2)=5.d0
      y(3)=4.d0*y(1)/(3.d0*y(1)-8.d0)
      tend=0.1108d0
      rtol=1.d-3
      ymax0=1.d0
      hmax=tend-t
      h=1.d-3
      nrw=150
      niw=30
      ijob(1)=0
      ijob(2)=0
      ijob(3)=0
      ijob(4)=0
      ijob(5)=0
      ijob(6)=2
      ijob(7)=1
      call limex(n,nzc,nzv,lsc,fcn,jac,t,y,tend,rtol,ymax0,hmax,h,ijob,
     1           nrw,rw,niw,iw)
      stop
      end
c
      subroutine fcn (n,nzv,t,y,dy,bv,irv,icv)
      implicit real*8(a-h,o-z)
      dimension y(n),dy(n),bv(1),irv(1),icv(1)
      dy(1)=-0.8d0*y(1)+1.d1*y(2)-0.6d0*y(1)*y(3)
      dy(2)=-1.d1+1.6d0*y(3)/y(2)
      dy(3)=0.8d0*y(1)+1.6d0*y(3)-0.6d0*y(1)*y(3)
      bv(1)=1.d0/y(2)
      irv(1)=2
      icv(1)=2
      return
      end
c
      subroutine lsc (nzc,bc,irc,icc)
      implicit real*8(a-h,o-z)
      double precision bc(nzc)
      dimension irc(nzc),icc(nzc)
      bc(1)=1.d0
      irc(1)=1
      icc(1)=1
      return
      end
c
      subroutine jac (n,mb,t,y,yp,a)
      implicit real*8(a-h,o-z)
      double precision y(n),yp(n),a(mb,n)
      a(1,1)=-0.8d0-0.6d0*y(3)
      a(1,2)=1.d1
      a(1,3)=-0.6d0*y(1)
      a(2,1)=0.d0
      a(2,2)=(yp(2)-1.6d0*y(3))/y(2)**2
      a(2,3)=1.6d0/y(2)
      a(3,1)=0.8d0-0.6d0*y(3)
      a(3,2)=0.d0
      a(3,3)=1.6d0-0.6d0*y(1)
      return
      end
      subroutine jdum (n,mb,t,y,yp,a)
      implicit real*8(a-h,o-z)
      return
      end

