program test_thomas
      implicit double precision (a-h,o-z)
      parameter (kdim=40)
      dimension bef(kdim),aft(kdim),rhsol(kdim),diag(kdim)
c
      n=6
      diag(1)=1.d0
      aft(1)=0.d0
      rhsol(1)=10.d0
        do i=2,n-1
          diag(i)=-2.d0
          bef(i) = 1.d0
          aft(i) = 1.d0
          rhsol(i)=0.d0
        enddo
      diag(n)=1.d0
      bef(n)=0.d0
      rhsol(n)=10.d0*n
      call trdiag (kdim,n,bef,diag,aft,rhsol)
      do i=1,n
      write(*,*)i,rhsol(i)
      enddo
c

      stop
      end

c
c#################################################################c
      subroutine trdiag (kdim,n,bef,diag,aft,rhsol)
c#################################################################c
c     solution of a tridiagonal matrix
      implicit double precision (a-h,o-z)
      dimension bef(kdim),aft(kdim),rhsol(kdim),diag(kdim)
c
      eps=1.d-8
      do 10 i=2,n
      if(dabs(diag(i-1)).lt.eps) stop  'Divide by zero'
      r = bef(i)/diag(i-1)
      diag(i)  = diag(i)  - r* aft(i-1)
      rhsol(i) = rhsol(i) - r* rhsol(i-1)
  10  continue
c
      rhsol(n) = rhsol(n)/diag(n)
      do 20 i=n-1,1,-1
      rhsol(i) = (rhsol(i)-aft(i)*rhsol(i+1)) / diag(i)
  20  continue
c
      return
      end
c