subroutine spline(r,start,n,db,db2)
c  spline calculates the cubic spline interpolation of the density
c  together with subroutine splint.
c  the main ideas are "stolen" from "numerical recipes", cambridge
c  university press (1992).
c  input r     : position coordinate (array)
c  input start : index, r(start) = starting point of the spline
c  input n     : index, r(n) = r_max
c  input db    : spin up and down density. (array)
c  output db2  : second derivative of db (array)
      implicit logical (a-z)
      double precision r,db,db2,u,sig,p,qn,un
      integer spin,start,n,i
      dimension r(561),db(2,561),db2(2,561),u(561)
      do 100 spin = 1,2
         db2(spin,start) = 0.d0
         u(start) = 0.d0
         do 10 i = start+1,n-1
            sig = (r(i)-r(i-1))/(r(i+1)-r(i-1))
            p = sig*db2(spin,i-1)+2.d0
            db2(spin,i) = (sig-1.d0)/p
            u(i) = (6.d0*((db(spin,i+1)-db(spin,i))/(r(i+1)-r(i))-
     j           (db(spin,i)-db(spin,i-1))/(r(i)-r(i-1)))/
     j           (r(i+1)-r(i-1))-sig*u(i-1))/p
 10      continue
         qn = 0.d0
         un = 0.d0
         db2(spin,n) = (un-qn*u(n-1))/(qn*db2(spin,n-1)+1.d0)
         do 20 i = n-1,start,-1
            db2(spin,i) = db2(spin,i)*db2(spin,i+1)+u(i)
 20      continue
 100  continue
      return
      end