c-----------------------------------------------------------
c     Integrates static equations of motion for spherically 
c     symmetric, general-relativistic boson star.
c
c     See class notes and Colpi et al, Phys Rev Lett,
c     vol 57, 2485--2488 for more details
c-----------------------------------------------------------
c     neq = 4
c-----------------------------------------------------------

      subroutine fcn(neq,x,y,yprime)
         implicit    none

         include    'fcn.inc'

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

         real*8      A,     M,     B,    sigma,    xi

         M     = y(1)
         B     = y(2)
         sigma = y(3)
         xi    = y(4)

         if( x .eq. 0.0d0 ) then
            A = 1.0d0 
            yprime(1) = 0.0d0
            yprime(2) = 0.0d0
            yprime(3) = 0.0d0

c---        The following expression previously had a sign error, but LSODA
c---        is good enough to ensure that results are apparently quite 
c---        insensitive to the gaffe.

            yprime(4) = - ( (Omegasq / B - 1.0d0) * sigma -
     &                     Lambda * sigma**3 ) / 3.0d0


         else
            A = 1.0d0 / (1.0d0 - 2.0d0 * M / x)
            yprime(1) = x**2 * (
     &          0.5d0 * (Omegasq / B + 1.0d0) * sigma**2 +
     &          0.25d0 * Lambda * sigma**4 + 
     &          0.5d0 * xi**2 / A
     &                          )
            yprime(2) = A * B * x * (
     &         (1.0d0 - 1.0d0 / A) / (x**2) + 
     &         (Omegasq / B - 1.0d0) * sigma**2 -
     &         0.5d0 * Lambda * sigma**4 +
     &         xi**2 / A
     &                               )
            yprime(3) = xi
            yprime(4) = xi * (
     &         -2.0d0 / x - 0.5d0 * yprime(2) / B + 
     &         A * (yprime(1) - M / x) / x 
     &                       )  - A * (
     &         (Omegasq / B - 1.0d0) * sigma - 
     &         Lambda * sigma**3      )
         end if

         return
      end 

c-----------------------------------------------------------
c     Implements Jacobian  (optional) ...
c-----------------------------------------------------------
      subroutine jac
         implicit    none

         include    'fcn.inc'

         return
      end