c  22. $B1"2rK!$G5a$a$kG.EAF3J}Dx<0$N(B3$B<!85%0%i%U(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex22
c
      implicit none
      integer :: tp,tp1,xp,k1,i,j,n
      real :: r,r1,r2,x,a(21,21),u(21,201)
c
      r =1.0
      r1=(2.0*(1.0+r))/r
      r2=(2.0*(1.0-r))/r
      k1 =4
      tp =51
      tp1=(tp-1)*k1+1
      xp =21
      n  =xp

      do i=1,xp
         x=0.05*float(i-1)
         if (i <= 11) then
           u(i,1)=2.0*x
         else
           u(i,1)=2.0*(1.0-x)
         endif
      enddo
      do j=1,tp1
         u( 1,j)=0.0
         u(21,j)=0.0
      enddo

      do i=1,n
         a(i,i)=r1
      enddo

      do i=2,n-1
         a(i,i+1)=-1.0
      enddo

      do i=2,n-1
         a(i,i-1)=-1.0
      enddo

         a(n,n-1)=-2.0      

      call inkaihou22(a,u,n,xp,tp1,r2,k1)

      stop 
      end
c
c-----$B1"2rK!$r2r$/(BLU$BJ,2rK!$rMQ$$$?%5%V%k!<%A%s(B
      subroutine inkaihou22(a,u,n,xp,tp,r2,k1)     
      implicit none
      integer :: i,j,k,n,m,l,maxi,xp,tp,time,k1
      real :: maxa,r2,xx,tt
      real :: a(n,n+1),aa(n,n),bb(n),ll(n,n),
     -        uu(n,n),da(n),db,ma(3),x(n),y(n),
     -        u(xp,tp),ub(n)
      
      do time=1,tp

      ub(1)= r2*u(1,time)+u(2,time)
      ub(n)=2.0*u(n-1,time)+r2*u(n,time)

      do i=2,n-1
         ub(i)=u(i-1,time)+r2*u(i,time)+u(i+1,time)
      enddo

      do i=1,n
       do j=1,n
        aa(i,j)=a(i,j)
        bb(i)=ub(i)
       enddo
      enddo

      do k=1,1
          maxi=k
          do i=k,n
             if ( abs(a(i,k)) > abs(a(maxi,k)) ) then
                maxi=i
             endif
          enddo
          do j=k,n
             da(j)     =aa(maxi,j)
             aa(maxi,j)=aa(k,j)
             aa(k,j)   =da(j)
             db        =bb(maxi)
             bb(maxi)  =bb(k)
             bb(k)     =db
          enddo
      enddo

      do k=1,n
         do j=1,n
          if (j <= k) then
            ll(k,j)=aa(k,j)
             do m=1,j-1
               ll(k,j)=ll(k,j)-ll(k,m)*uu(m,j)
             enddo
          elseif (j > k) then
            uu(k,j)=aa(k,j)/ll(k,k)
             do m=1,k-1
                uu(k,j)=uu(k,j)-(ll(k,m)*uu(m,j))/ll(k,k)
             enddo
          endif
         enddo
      enddo
      do i=1,n
       do j=i,n
          aa(i,j)=ll(i,j)+uu(i,j)
       enddo
      enddo
c      write(6,*)''
c      write(6,*)'LU$BJ,2r(B:'
c      write(6,*)''
c      write(6,*)'Ly=b:'
      do i=1,n
         y(i)=bb(i)/ll(i,i)
         do j=1,i-1
            y(i)=y(i)-(ll(i,j)*y(j))/ll(i,i)
         enddo
      enddo
c      write(6,'(a4,f7.4)')'y1=',y(1)
c      write(6,'(a4,f7.4)')'y2=',y(2)
c      if (n==3) then
c      write(6,'(a4,f7.4)')'y3=',y(3)
c      endif

c      write(6,*)''
c      write(6,*)'Ux=y:'
      do i=n,1,-1
         x(i)=y(i)
         do j=n,i+1,-1
            x(i)=x(i)-uu(i,j)*x(j)
         enddo
      enddo

      do i=2,n-1
         u(i,time+1)=x(i)
      enddo

      enddo    !$B$+$J$j>e$N!V(Btime$B!W$KBP1~(B

      do j=1,tp,k1
         do i=1,xp
            tt=(0.01/float(k1))*float(j-1)
            xx=0.05*float(i-1)
            write(10,*) tt,xx,u(i,j)
         enddo
         write(10,*) ''
      enddo


c      write(6,'(a4,f7.4)')'x1=',x(1)
c      write(6,'(a4,f7.4)')'x2=',x(2)
c      if (n==3) then
c      write(6,'(a4,f7.4)')'x3=',x(3)
c      endif

c      write(6,*)''
c      write(6,*)'$B2r(B'
c      write(6,'(a3,f7.4)')'x=',x(1)
c      write(6,'(a3,f7.4)')'y=',x(2)
c      if (n==3) then
c      write(6,'(a3,f7.4)')'z=',x(3)
c      endif
c      go to 30
c   10 write(6,*)''
c      write(6,*)'$B2r$OITDj(B'  
c      go to 30
c   20 write(6,*)'$B2r$OITG=(B'
c   30 continue  
      return
  100 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,1x,f7.4,1x,f7.4) 
  200 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,
     -       1x,a3,1x,f7.4) 
  300 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,1x,f7.4,
     -       1x,a3,1x,f7.4) 
      end

