c-----------------------------------------------------------------------
c     these are the subroutines for the
c
c        VANDERPOL 6
c        ODE of dimension 6
c
c-----------------------------------------------------------------------
      subroutine prob(problm,neqn,tbegin,tend,ijac,mljac,mujac)
      character*(*) problm
      integer neqn,ijac,mljac,mujac
      double precision tbegin,tend,x
      common /final/ x

      problm = 'vdpol6'
      neqn   = 6
      tbegin = 0d0
      tend   = 20.0d0
      ijac   = 1
      mljac  = neqn
      mujac  = neqn

      x = tend
      return
      end
c-----------------------------------------------------------------------
      subroutine init(neqn,y)
      integer neqn
      double precision y(neqn)

c     integer i
        do  i=1,6
           y(i)=1.0D+0
        enddo
      return
      end
c-----------------------------------------------------------------------
        SUBROUTINE feval(N,X,Y,DF)
C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),DF(N)
        df(1)=-10.0D+0*y(1)+1000.0D+0*y(2)
        df(2)=-1000.0D+0*y(1)-10.0D+0*y(2)
        df(3)=-4.0D+0*y(3)
        df(4)=-y(4)
        df(5)=-0.5D+0*y(5)
        df(6)=-0.1D+0*y(6)
        RETURN
        END

c-----------------------------------------------------------------------
      subroutine jeval(N,X,Y,DFY,MEBND)
C --- JACOBIAN OF VANDERPOL EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),DFY(MEBND,N)
        DFY(1,1)=-10.0D0
        DFY(1,2)=1000.D0
        DFY(2,1)=-1000.0D+0
        DFY(2,2)=-10.0D+0
        dfy(3,3)=-4.0D+0
        dfy(4,4)=-1.0D+0
        dfy(5,5)=-0.5D+0
        dfy(6,6)=-0.1D+0
        RETURN
        END

c-----------------------------------------------------------------------
      subroutine solut(neqn,true)
      integer neqn
      double precision true(neqn), x
      common /final/ x
c
c
         true(1)=dexp(-10.0D+0*x)*(cos(1000.0D+0*x)+sin(1000.0D+0*x))
         true(2)=dexp(-10.0D+0*x)*(cos(1000.0D+0*x)-sin(1000.0D+0*x))
         true(3)=dexp(-4.0D+0*x)
         true(4)=dexp(-x)
         true(5)=dexp(-0.5D+0*x)
         true(6)=dexp(-0.1D+0*x)

      return
      end