c-----------------------------------------------------------------------
c     these are the subroutines for the
c
c        Robertson
c        ODE of dimension 3
c
c-----------------------------------------------------------------------

      subroutine prob(problm,neqn,tbegin,tend,ijac,mljac,mujac)
      character*(*) problm
      integer neqn,ijac,mljac,mujac
      double precision tbegin,tend

      problm = 'rober'
      neqn   = 3
      tbegin = 0d0
      tend   = 4d6
      ijac   = 0
      mljac  = neqn
      mujac  = neqn

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

      y(1) = 1d0
      y(2) = 0d0
      y(3) = 0d0

      return
      end
c-----------------------------------------------------------------------
      subroutine feval(neqn,t,y,dy)
      integer neqn
      double precision t,y(neqn),dy(neqn)
      dy(1) = -.04d0*y(1) + 1.d4*y(2)*y(3)
      dy(3) = 3.d7*y(2)*y(2)
      dy(2) = -dy(1) - dy(3)
      return
      end
c-----------------------------------------------------------------------
      subroutine jeval(neqn,t,y,jac,ldim)
      integer neqn,ldim
      double precision t,y(neqn),jac(ldim,neqn)

      integer i,j

      do 20 j=1,neqn
         do 10 i=1,neqn
            jac(i,j)=0d0
   10    continue
   20 continue

      jac(1,1) = -.04d0
      jac(1,2) = 1.d4*y(3)
      jac(1,3) = 1.d4*y(2)
      jac(2,1) = .04d0
      jac(2,3) = -jac(1,3)
      jac(3,1) = 0.0
      jac(3,2) = 6.d7*y(2)
      jac(3,3) = 0.0
      jac(2,2) = -jac(1,2) - jac(3,2)
      return
      end
c-----------------------------------------------------------------------
      subroutine solut(neqn,y)
      integer neqn
      double precision y(neqn)
      y(1) =  .51680960148711D-03
      y(2) =  .20682944912031D-08
      y(3) =  .99948318833021D+00

      return
      end