c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 9.2
c
c     Example of solving a P.D.E using implicit method
c
c
c     file:  exs92.f
c
      parameter (n=10)
      dimension v(n),a(n),d(n),c(n),error(n)
      data M,hk/200,0.005125/  
c
      print *
      print *,' Implicit method example'
      print *,' Section 9.2, Kincaid-Cheney'
      print *
      print 7
c
      pi = 4.0*atan(1.0)
      h = 1.0/real(n+1)
      s = hk/h**2 
c
      do 2 i=1,n
         v(i) = sin(pi*real(i)*h) 
         d(i) = 1.0 + 2.0*s
 2    continue
c
      t = 0.0
      print 8,0,t,v(0),v(1),v(10),v(11)
      print *
c
      do 3 i=1,n-1
         c(i) = -s
         a(i) = -s
 3    continue 
c
      do 6 j=1,M
         call tri(n,a,d,c,v,v)
         t = real(j)*hk
         if (0 .eq. mod(j,10)) then 
            print 8,j,t,v(0),v(1),v(10),v(11)
            print *
         endif
c
         do 5 i=1,n
            x = real(i)*h
            error(i) = abs(sin(pi*x)/exp(pi*pi*t) - v(i))
 5       continue
c   
         if (0 .eq. mod(j,10)) then
            print *,' norm of error =', unorm(n,error)
            print *
         endif
 6    continue
c
 7    format (3x,'j',8x,'t',13x,'v(0)',11x,'v(1)',6x,'.....',
     +        5x,'v(10)',10x,'v(11)')
 8    format (1x,i3,2x,3(e13.6,2x),5x,2(e13.6,2x))
      stop
      end
c
      subroutine tri(n,a,d,c,b,x)     
c
c     solve a tridiagonal system
c
      dimension a(n),d(n),c(n),b(n),x(n) 
c
      do 2 i = 2,n
         xmult = a(i-1)/d(i-1) 
         d(i) = d(i) - xmult*c(i-1)    
         b(i) = b(i) - xmult*b(i-1)    
 2    continue    
c
      x(n) = b(n)/d(n)      
      do 3 i = n-1,1,-1     
         x(i) = (b(i) - c(i)*x(i+1))/d(i)
 3    continue    
c
      return      
      end 
c
      function unorm(n,u)
c
c     computes infinity norm of n component vector v
c
      dimension u(n)
      real max
c
      temp = 0.0
      do 2 i=1,n
         temp = max(temp,abs(u(i)))
 2    continue
      unorm = temp
c
      return
      end