program sturm
      real   nterm, nx
      parameter (nterm = 5, nx = 11)
c
c  sturm computes the solution of the Sturm-Liouville problem,
c  p = q = 1, f = -sig(al*sin(l-0.5)*pi*x, using the finite element
c  method on a uniform grid with either linear or quadrtic 
c  interpolation.
c  the banded system is solved with banfac, bansol 
c  
      real      y(nx),yex(nx),x(nx),b(5,nx),g(nx),f(nx),
     +          a(nterm),yd(nx)
c
      open(1,file='sturm.dat')      
      open(6,file='sturm.out')
      read(1,1) int,ipr
      do 60 i = 1,nterm
   60 read(1,2) a(i)
    1 format(i8)
    2 format(e10.3)
c
      if (int .eq. 1) write(6,3)
      if (int .eq. 2) write(6,4)
      write(6,5) nx,a
    3 format(' Sturm-Liouville problem, fem: linear interpolation')
    4 format(' Sturm-Liouville problem, fem: quadratic interpolation')
    5 format(' nx=', i3,'  a=', 5e10.3,//)
c      
      pi = 3.1415927
      nxp = nx-1
      dx = 1./(nx-1)
      dxs = dx*dx
c
c   set grid, exact solution and rhs(f)
c      
      do 6 i = 1,nx
       ai = i-1
       x(i) = ai*dx
       f(i) = 0.
       yex(i) = 0.
       do 6 l = 1,nterm
        al = (l-0.5)*pi
        dum = a(l)*sin(al*x(i))
        f(i) = f(i) - dum
        yex(i) = yex(i) - dum/(1.-al*al)
    6 continue
c
c   set up arrays, b and g
c      
      if(int .eq. 1) then
       do 7 i = 2,nx,int
        im = i-1
        b(1,im) = 0
        b(2,im) = 1./dxs + 1./6.
        b(3,im) = -2./dxs + 2./3.
        b(4,im) = b(2,im)
        b(5,im) = 0.
        g(im)   = ( f(i-1) + 4.*f(i) + f(i+1) )/6.
    7  continue
      else
       do 8 i = 2,nx,int
        im = i-1
        b(1,im) = 0
        b(2,im) = 4.*( 1./dxs + 0.1 )/3.
        b(3,im) = 4.*( -2./dxs + 0.8 )/3.
        b(4,im) = b(2,im)
        b(5,im) = 0.
        g(im)   = 0.4*( f(i-1) + 8.*f(i) + f(i+1) )/3.
        b(1,i) = 2.*( -.25/dxs - 0.1 )/3.
        b(2,i) = 2.*( 2./dxs + 0.2 )/3.
        b(3,i) = 2.*( -3.5/dxs + 0.8 )/3.
        b(4,i) = b(2,i)
        b(5,i) = b(1,i)
        g(i)   = -0.1*( f(i-1)+f(i+3) ) + 0.2*( f(i)+f(i+2) ) 
     +          + 0.8*f(i+1)
        g(i)   = 2.*g(i)/3.
    8  continue
      endif
c
c   adjust for boundary conditions
c      
      b(2,1) = 0.
      b(1,2) = 0.
      b(3,nxp) = 0.5*b(3,nxp)
      b(4,nxp) = 0.
      b(5,nxp) = 0.
      if (int .eq. 1) g(nxp) = ( f(nxp)+2.*f(nx) )/6.
      if (int .eq. 2) 
     +     g(nxp) = 2.*( -0.1*f(nxp-1)+0.2*f(nxp)+0.4*f(nx) )/3.
c
c   solve banded system of equations
c      
      call banfac(b,nxp,int)
      
      call bansol(g,yd,b,nxp,int)
c
c   compute rms error
c      
      sum = 0.
      do 10 i = 2,nx
       y(i) = yd(i-1)
       dy = y(i) - yex(i)
       sum = sum + dy*dy
       if (ipr .eq. 0) goto 10
       write(6,9) i, x(i), y(i), yex(i), dy
    9  format('  i=',i3,'  x=',f7.5,'  y=',f7.5,'  yex+',f7.5,
     +        '  dy=',f7.5)
   10 continue
      rms =sqrt(sum/(nx-1))
      write(6,11) rms,nx
   11 format(//,' rms=',e10.3,'   nx=',i3)
c      
      stop
      end 
c ------------ : ------------
      subroutine banfac(b,n,int)
      real   nterm, nx
      parameter (nterm = 5, nx = 11)
c
c  factorizes band matrix arising from linear or quadratic elements
c  into l.u
c  
      real     b(nterm,nx)
      
c
c  int = 1, linear elements = tridiagonal system
c  int = 2, qudratic elements = pentadiagonal system
c  assumes that first equation formed at midside node
c  
      if (int .eq. 1) then
       np = n-1
       do 1 j = 1,np
        jp = j+1
        b(2,jp) = b(2,jp)/b(3,j)
        b(3,jp) = b(3,jp) - b(2,jp)*b(4,j)
    1  continue
       return
      else
       nh = n/2
       do 5 i = 1,2
        js = 3-i
        do 4 j = js,nh
         ja = 2*(j-1)
c  i = 1, first pass, reduction to tridiagonal
c  i = 2, second pass, reduction to upper triangular
         if (i .eq. 1) then
          jb = ja+2
          b(1,jb) = b(1,jb)/b(2,jb-1)
          b(2,jb) = b(2,jb) - b(1,jb)*b(3,jb-1)
          b(3,jb) = b(3,jb) - b(1,jb)*b(4,jb-1)
         else
          jb = ja+3
          b(2,jb-1) = b(2,jb-1)/b(3,jb-2)
          b(3,jb-1) = b(3,jb-1) - b(2,jb-1)*b(4,jb-2)
          if (jb .le. n) then
           b(2,jb) = b(2,jb)/b(3,jb-1)
           b(3,jb) = b(3,jb) - b(2,jb)*b(4,jb-1)
           b(4,jb) = b(4,jb) - b(2,jb)*b(5,jb-1)
          endif
         endif
    4   continue
    5  continue
      endif
      return
      end
c ------------ : ------------
      subroutine bansol(r,x,b,n,int)
      real   nterm, nx
      parameter (nterm = 5, nx = 11)
c
c  uses l.u factorization to solve for x, given r
c  
      real     r(nx),x(nx),b(5,nx)

c
c  int = 1, tridiagonal system
c  int = 2, pentadiagonal system
c   assumes that first equation formed at midside node
c   ibc=0 if last equation formed at midside node
c   ibc=1 if last equation formed at corner node
c  
      if (int .eq. 1) then
       np = n-1
       do 1 j = 1,np
        jp = j+1
        r(jp) = r(jp) - b(2,jp)*r(j)
    1  continue
       x(n)=r(n)/b(3,n)
       do 2 j = 1,np
        ja = n - j
        x(ja) = ( r(ja) - b(4,ja)*x(ja+1) )/b(3,ja)
    2  continue
       return
      else
    3  ibc = 1
       nh = n/2
       if (2*nh .eq. n) r(n+1) = 0.
       if (2*nh .eq. n) b(2,n+1) = 0.
       do 6 i = 1,2
        do 5 j = 1,nh
         ja = 2*j
         do 4 k = 1,i
          jb = ja -1 + k
          r(jb) = r(jb) - b(i,jb)*r(jb-1)
    4    continue
    5   continue
    6  continue
       nen = nh -ibc
       x(n) = r(n)/b(3,n)
       if ( ibc .eq. 1) x(n-1) = ( r(n-1) - b(4,n-1)*x(n) )/b(3,n-1)
       do 7 j = 1,nen
        ja = n -2*j + 1 - ibc
        x(ja) = ( r(ja) - b(4,ja)*x(ja+1) - b(5,ja)*x(ja+2) )/b(3,ja)
        x(ja-1) = ( r(ja-1) - b(4,ja-1)*x(ja) )/b(3,ja-1)
    7  continue
      endif
      return
      end
c ------------ : ------------