subroutine thomas(a,b,c,r,n)
c
c  Tri-diagonal matrix solver
c
c
c   Be sure of your precision
c     implicit double precision (a-h,o-z)
c
      dimension a(n),b(n),c(n),r(n)
      c(1) = c(1)/b(1)
      r(1) = r(1)/b(1)
      do 40 i = 2,n
       im = i-1
       c(i) = c(i)/(b(i)-a(i)*c(im))
       r(i) = (r(i)-a(i)*r(im))/(b(i)-a(i)*c(im))
  40  continue
      do 50 i = 1,n-1
       j = n-i
       jp = j+1
       r(j) = r(j) - c(j)*r(jp)
   50 continue
      return
      end