program thomas
      implicit double precision(a-h,o-z)
      parameter(NDIM=1001)
      dimension u(NDIM),phi(NDIM),phip(NDIM)
      data acc/1.0d-10/
      pi = acos(-1d0)

      write(6,1000) '  Enter z and n  : '
 1000 format(/a,$)
      read(5,*) z,n

*     xi = sqrt(9d0*pi**2/(128d0*z))
*    correction by G. Faussurier (7/4/04)
      xi = (9d0*pi**2/(128d0*z))**(1d0/3d0)
*

      xdphi = -(z-n)/z
 100  continue
      write(6,1000) '  Enter two initial guesses for R : '
      read(5,*) ra,rb
      a = ra/xi
      b	= rb/xi
      i = 0
 200  continue
         i = i + 1
         if(i.eq.1) then
            x0 = a
         else if(i.eq.2) then
            x0 = b
         else
            if (fa*fb.gt.0.0d0) then
               write(6,*) '  No Solution with these guesses !'
               write(6,1010) '  Ra = ',xi*a,'  phi(0,a)-1 = ',fa
 1010          format(a,f10.6,a,f16.8)
               write(6,1010) '  Rb = ',xi*b,'  phi(0,b)-1 = ',fb
               go to 100
            else
               x0 = 0.5*(a+b)
            endif
         endif
         u(1) = sqrt(x0)
         h = -u(1)/(NDIM-1)
         do 210 j= 1, NDIM
            u(j) = u(1) + (j - 1) * h
 210     continue
         phi(1) = 0
         phip(1) = 2*xdphi/u(1)
         call runge(u,phi,phip,h)
         write(6,1020) ' Initial value for R =',xi*x0
 1020    format(/a,f16.10)
         write(6,1030) '  phi(0) = ',phi(NDIM)
 1030    format(a,f14.10/)
         if(i.eq.1) fa = phi(NDIM) - 1.0d0
         if(i.eq.2) fb = phi(NDIM) - 1.0d0
         if(i.gt.2) then
             fc = phi(NDIM) - 1.0d0
             if (fc*fa.lt.0.0d0) then
                 b = x0
                 fb = fc
             else
                 a = x0
                 fa = fc
             endif
         endif
      if(abs(phi(NDIM) - 1.0d0).gt.acc) go to 200

*****  output file (1)
*
*   thomas.dat: r[i], phi[i], N[i], Zeff[i]
*
*   where V(r) == -Z/r + Vel(r) = -Zeff(r)/r 
*
*****       
      rmax = xi*u(1)**2
      open(unit=1,file='thomas.dat',form='formatted',status='unknown')
      
      rewind(1)
	write(1,1090)
 1090 format(7x,'r',10x,'phi',10x,'N',10x,'Zeff')
     
      do 300 i = NDIM,1,-2
          r = xi*u(i)**2
          ph = phi(i)
          en = z*(0.5d0*u(i)*phip(i)-phi(i)+1)
          zeff = z*ph + (z-n)*r/rmax        
          write(1,1100) r,ph,en,zeff           
 1100     format(4f12.5)
 
 300  continue 
      rewind(1)
      stop
      end 

      double precision function f(x,y,yp)
      implicit double precision(a-h,o-z)
      parameter(NDIM=1001)

      if (x.eq.0.0d0) then
          f = 0.0d0
      else
          f = 4.0 * x * sqrt(abs(y*y*y)) + yp / x
      endif
      return
      end

      subroutine runge(x,y,yp,h)
      implicit double precision(a-h,o-z)
      double precision k1,k2,k3,k4
******************************************************************************
*    Fourth-Order Runge-Kutta Method
*    Abramowitz and Stegun 
*    Handbook of Math Functions (1964 edition)
*    NBS AMS 55
*    Chapter 25, Equation 25.5.20
******************************************************************************
      parameter(NDIM=1001)
      dimension x(NDIM),y(NDIM),yp(NDIM)
      data six/6.0d0/,two/2.0d0/,half/0.5d0/,eighth/0.125d0/
      do 100 i = 1,NDIM-2
          u = x(i)
          v = y(i)
          w = yp(i)
          w = yp(i)
          k1 = h * f(u,v,w)
          u = x(i) + half * h
          v = y(i) + h * (half * yp(i) + eighth * k1)
          w = yp(i) + half * k1
          k2 = h * f(u,v,w)
          u = x(i) + half * h
          v = y(i) + h * (half * yp(i) + eighth * k2)
          w = yp(i) + half * k2
          k3 = h * f(u,v,w)
          u = x(i) + h
          v = y(i) + h * ( yp(i) + half * k3 )
          w = yp(i) +k3
          k4 = h * f(u,v,w)
          y(i+1) = y(i) + h * ( yp(i) + ( k1 + k2 + k3 )/six )
          yp(i+1) = yp(i) + (k1 + two * k2 + two * k3 + k4)/six
 100  continue
*
*        Special treatnent for x = 0.0
*        Adams-Bashforth Fourth-Order 
*        Chapter 25, Equation 25.5.4
*
         i = NDIM-1
         yp(i+1) =  yp(i) +
     &             h * ( 55.0d0 * f(x(i),y(i),yp(i))
     &             -59.0d0 * f(x(i-1),y(i-1),yp(i-1))
     &             +37.0d0 * f(x(i-2),y(i-2),yp(i-2))
     &             - 9.0d0 * f(x(i-3),y(i-3),yp(i-3)) ) /24.0d0
         y(i+1)  =  y(i) +
     &                   h * ( 55.0d0 * yp(i)
     &                   -59.0d0 * yp(i-1)
     &                   +37.0d0 * yp(i-2)
     &                   - 9.0d0 * yp(i-3)  ) /24.0d0

      return
      end