c
c vanderpol_rhs.f
c
c Vector field functions for the vector field 'vanderpol'
c These functions are to be used with the Fortran ODE solver RADAU5.
c
c This file was generated by the program VFGEN (Version:2.4.0)
c Generated on 10-Jul-2008 at 17:47
c
      subroutine vanderpol_rhs(n_,t_,y_,f_,rpar_,ipar_)
      implicit none
      integer n_, ipar_
      double precision t_, y_, f_, rpar_
      dimension y_(2), f_(2), rpar_(1)
      double precision epsilon
      double precision x, y

c     --- Parameters ---
      epsilon    = rpar_(1)
c     --- State variables ---
      x          = y_(1)
      y          = y_(2)
c     --- The vector field ---
      f_(1) = 1.0D0/epsilon*( y+x+-3.3333333333333331D-01*x**3.0D0)
      f_(2) = -x

      return
      end

c
c Jacobian of the vector field 'vanderpol'
c
      subroutine vanderpol_jac(n_,t_,y_,dfy_,ldfy_,
     &                             rpar_,ipar_)
      implicit none
      integer n_, ldfy_, ipar_
      double precision t_, y_, dfy_, rpar_
      dimension y_(2), dfy_(ldfy_,2), rpar_(1)
      double precision epsilon
      double precision x, y

c     --- Parameters ---
      epsilon    = rpar_(1)
c     --- State variables ---
      x          = y_(1)
      y          = y_(2)
c     --- Jacobian ---
      dfy_(1, 1) = 1.0D0/epsilon*(-x**2.0D0+1.0D0)
      dfy_(1, 2) = epsilon**-1.0D0
      dfy_(2, 1) = -1.0D0
      dfy_(2, 2) = 0.0D0

      return
      end

c
c vanderpol_out -- output subroutine
c
      subroutine vanderpol_out(nr,told,t,y,cont,lrc,n,
     &                  rpar,ipar,irtrn)
      implicit none
      integer nr, lrc, n, ipar, irtrn
      double precision told, t, y, cont, rpar
      integer i
      dimension y(2), rpar(1)

      write (6,99) t, (y(i), i = 1,2)
99    format(1x,f10.5,2E18.10)
      return
      end