      subroutine plfe1d(n,u,v,ds,us,dm,um,x,a,rm,init)
*******************************************
*  Generic piece-wise linear FE solver for 1-D equation:
*   rm *u -(a(x)u_x)_x =f, u(x(0))=u(x(n))=0.
* Input: 
*     u(i),i=1,n-1: RHS, unchanged at output
*     x(i),i=0,n:  The collocation points in ascending order
*     a(i),i=0,n: The variable coefficients (a(i) >= 0)
*     init.eq.1: Initialize the mass and stiffness matrices
* Output:
*     v(i), i=1,n-1: The solution
*     dm,um: The diagonal and upper diagonal of the Mass matrix
*     ds, us: The diagonal and upper diagonal of (rm*Mass.+Stiff.) 
*             after calling dpttrs.f. Can be repeated used in dpttrs.f!
*
* Rmk. In the integration, \int_{x_j}^{x_{j+1}} a(x) f(x) dx is replaced 
*      by  .5d0*(a(j)+a(j+1))*\int_{x_j}^{x_{j+1}} f(x) dx
*******************************************
      implicit double precision (a-h,o-z)
      dimension u(n-1),v(n-1),x(0:n),ds(n-1),us(n-1),dm(n-1),um(n-1)
      dimension a(0:n)

      if (x(0).gt.x(1)) then
         print*,'Error: collocation points should be'
         print*,'       in ascending order  in plfe1d.f !'
         stop
      endif
      if (init.eq.1) then
**** setup the matrix 

**** diagonal elements 
      do j=1,n-1
         t1=.5d0*(a(j)+a(j-1))
         t2=.5d0*(a(j)+a(j+1))
         dm(j)=t1*(x(j)-x(j-1))/3+t2*(x(j+1)-x(j))/3
         ds(j)=t1/(x(j)-x(j-1))+t2/(x(j+1)-x(j))
         ds(j)=ds(j)+rm*dm(j)
      enddo

***** Off-diagonal elements
      do j=1,n-2
         t2=.5d0*(a(j)+a(j+1))
         um(j)=t2*(x(j+1)-x(j))/6
         us(j)=-t2/(x(j+1)-x(j))
         us(j)=us(j)+rm*um(j)
      enddo

      call dpttrf(n-1,ds,us,info)
      endif

***** multiply the righthand-side by the mass matrix
      do j=2,n-2
         v(j)=u(j)*dm(j)+u(j-1)*um(j-1)+u(j+1)*um(j)
      enddo
      v(1)=u(1)*dm(1)+u(2)*um(1)
      v(n-1)=u(n-1)*dm(n-1)+u(n-2)*um(n-2)
      call dpttrs(n-1,1,ds,us,v,n-1,info)

      return
      end
