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.6
c
c     Example of use of method of characteristics   
c
c 
c     file: ex3s96.f
c
      parameter(n=8)
      dimension x(n),y(n),u(n),p(n),q(n)
c
      print *
      print *,' Method of Characteristics'
      print *,' Section 9.6, Kincaid-Cheney'
      print *
      print 5
c
      h = 1.0/(n-1)
      do 2 j=1,n 
         x(j) = (j-1)*h
         y(j) = 0.0
         u(j) = f(x(j))
         p(j) = df(x(j))
         q(j) = g(x(j))
         print 6,x(j),y(j),p(j),q(j),u(j)
 2    continue
      print *
c
      do 4 i=1,n
         y(i) = i*h
         do 3 j=1,n-i
            x(j) = x(j) + h/2
            pp = (p(j) + p(j+1))/2. + (1.0 - h/8.)*(q(j+1) - q(j))
            qq = (p(j+1)-p(j) + (2.- h/4.)*(q(j+1) + q(j)))/(4.+ h/2.)
            u(j) = u(j) + h*(pp+p(j))/4. + h*(qq+q(j))/2.
            p(j) = pp
            q(j) = qq
            print 6,x(j),y(j),p(j),q(j),u(j)
 3       continue
         print *
 4    continue
c
 5    format (8x,'x',14x,'y',14x,'p',14x,'q',14x,'u')
 6    format (2x,5(e13.6,2x))
      stop
      end
c
      real function f(x) 
      f = exp(x/4.)
      return
      end
c
      real function g(x)
      g = (-1 + sqrt(2.))*exp(x/4.)/8.
      return
      end
c
      real function df(x)
      df = exp(x/4.0)/4.
      return
      end
c
      real function tu(x,y)
      tu = exp(x/4.)*exp((-1+sqrt(2.))*y/8.)
      return
      end