c===========================================================
c     Driver routine which integrates ODEs defining 
c     initial data for 1d, massive, complex scalar field in
c     polar-areal coordinates.
c 
c===========================================================
      subroutine fcn(neq,r,y,yprime)
         implicit    none

c-----------------------------------------------------------
c     Debug auxiliary variables 
c-----------------------------------------------------------
      character*(3)   cdnm
      parameter       ( cdnm = 'fcn' )

      logical         ltracet
      parameter       ( ltracet = .true. )

      logical         ltracef
      parameter       ( ltracef = .false. )

c-----------------------------------------------------------
         include    'fcn.inc'

         integer     neq
         real*8      r,     y(neq),    yprime(neq)

         real*8      a,     alpha
         real*8      u,     dudphi2

c-----------------------------------------------------------

       if(neq.eq.4)then
         a = y(1)
         alpha = y(2)
         phj = y(3)
         ppj = y(4)
       else
         a = y(1)
         alpha = y(2)
       endif

         call pot(neq,y,u,dudphi2)

       if ( ltracef ) then
        write (0,*) cdnm, ': check if the initial values are correct'
        write (0,*) cdnm, ': r = ', r 
        write (0,*) cdnm, ': neq = ', neq 
        write (0,*) cdnm, ': a = ', a
        write (0,*) cdnm, ': alpha = ', alpha 
        write (0,*) cdnm, ': phj = ', phj 
        write (0,*) cdnm, ': ppj = ', ppj 
        write (0,*) 
       endif

         if ( r .eq. 0.0d0) then 
         yprime(1) = 0.0d0 
         yprime(2) = 0.0d0 
          if(neq.eq.4)then
           yprime(3) = ppj

           yprime(4) =  - ( w**2/alpha**2 - dudphi2 ) * phj * a**2
          endif

         else
         yprime(1) = 0.5d0 *( a*(1.0d0-a**2)/r + 4.0d0*pie*r*a
     &                    *(a**2 * u +phj**2 * a**2 * w**2/alpha**2
     &                    +ppj**2)) 

         yprime(2) = alpha * 0.5d0 *( (a**2 - 1.0d0)/r + 4.0d0*pie*r
     &                    *(phj**2 * a**2 * w**2/alpha**2 - a**2 * u 
     &                    +ppj**2)) 

         if(neq.eq.4)then
         yprime(3) = ppj

         yprime(4) = - (1.0d0 +a**2 - 4.0d0*pie*r**2 * a**2 * u)*ppj/r
     &                       - ( w**2/alpha**2 - dudphi2 ) * phj * a**2
         endif

         endif

       if ( ltracef ) then
        write (0,*) cdnm, ': yprime = ', yprime
       endif

         return
      end 

c===========================================================
c     Dummy Jacobian routine.
c===========================================================
      subroutine jac
         implicit    none

         include    'fcn.inc'

         return
      end