c  22(2). $B1"2rK!$G5a$a$?CM$H2r@O2r$H$N8m:9(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex22a
c
      implicit none
      integer :: tp,tp1,xp,k1,i,j,n
      real :: r,r1,r2,x,t,a(21,21),u(21,201),gosa,
     -        ui(10,10),ut(10,10)
c
      write(6,*)'22(2). $B1"2rK!$G5a$a$?CM$H2r@O2r$H$N8m:9(B'
      write(6,*)''
      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 kaisekikai22a(ut)

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

      do j=21,201,20
         do i=2,11
            ui(i-1,(j-1)/20)=u(i,j)
         enddo
      enddo

      write(6,200)
      do j=1,10
         do i=1,10
            t=j*0.05 
            x=i*0.05
            gosa=abs(ui(i,j)-ut(i,j))
            write(6,100) t,x,ui(i,j),ut(i,j),gosa
         enddo
         write(6,*)''
      enddo

      stop 
  100 format(f5.2,1x,f5.2,e18.9,1x,e18.9,1x,e18.9) 
  200 format(4x,'t',5x,'x',3x,'$B1"2rK!$G$N2r(B',
     -       7x,'$B2r@O2r(B',13x,'$B@dBP8m:9(B') 
      end
c
c-----$B2r@O2r$N%5%V%k!<%A%s(B
      subroutine kaisekikai22a(ut)
      implicit none
      integer :: n0,m0,l0,ll,mm
      real :: ut(10,10),x,t,n,m,l,pi,tt,xx
c
      pi=3.141592653589
      mm=0
      do m0=1,10
        ll=0
        mm=mm+1
        m=float(m0)
        t=0.05*m
        do l0=1,10
          ll=ll+1
          l=float(l0)
          x=0.05*l
          ut(ll,mm)=0.0
          do n0=1,100
             n=float(n0)
             ut(ll,mm)=ut(ll,mm)+(8.0/(pi**2))*(1.0/(n**2))*
     -                 sin((n*pi)/2.0)*sin(n*pi*x)*
     -                 exp(-(n**2)*(pi**2)*t)      
          enddo
c        tt=0.01*m
c        xx=0.05*l
c        write(6,*) tt,xx,ut
c        write(6,*) 'ut(',ll,',',mm,')', ut(ll,mm)
        enddo
c      write(6,*)
      enddo      
      
      return
      end      

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
cc      write(6,*)''
cc      write(6,*)'LU$BJ,2r(B:'
cc      write(6,*)''
cc      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
cc      write(6,'(a4,f7.4)')'y1=',y(1)
cc      write(6,'(a4,f7.4)')'y2=',y(2)
cc      if (n==3) then
cc      write(6,'(a4,f7.4)')'y3=',y(3)
cc      endif

cc      write(6,*)''
cc      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

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


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

cc      write(6,*)''
cc      write(6,*)'$B2r(B'
cc      write(6,'(a3,f7.4)')'x=',x(1)
cc      write(6,'(a3,f7.4)')'y=',x(2)
cc      if (n==3) then
cc      write(6,'(a3,f7.4)')'z=',x(3)
cc      endif
cc      go to 30
cc   10 write(6,*)''
cc      write(6,*)'$B2r$OITDj(B'  
cc      go to 30
cc   20 write(6,*)'$B2r$OITG=(B'
cc   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

