c  21(2). $BM[2rK!$G5a$a$?CM$H2r@O2r$H$N8m:9(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex21a
c
      implicit none
      integer :: tp,tp1,tp2,xp,k1,k2,i,j
      real :: r1,r2,uy(10,10),ut(10,10),gosa,t,x
c
      write(6,*)'21(2). $BM[2rK!$G5a$a$?CM$H2r@O2r$H$N8m:9(B'
      write(6,*)''

      r1 =0.25
      r2 =0.5
      k1 =80
      k2 =40
      tp =11
      tp1=(tp-1)*k1+1
      tp2=(tp-1)*k2+1
      xp =21

      call kaisekikai21a(ut)

      write(6,*)'k=1/1600$B$G7W;;$7$?>l9g(B'
      call youkaihou21a(r1,k1,tp1,xp,uy)

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

      write(6,*)'k=1/800$B$G$G7W;;$7$?>l9g(B'
      call youkaihou21a(r2,k2,tp2,xp,uy)

      write(6,200)
      do j=1,10
         do i=1,10
            t=j*0.05 
            x=i*0.05
            gosa=abs(uy(i,j)-ut(i,j))
            write(6,100) t,x,uy(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,'$BM[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 kaisekikai21a(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-----$BM[2rK!$N%5%V%k!<%A%s(B
      subroutine youkaihou21a(r,k,tp,xp,uy)
      implicit none
      integer :: i,j,l,m,tp,xp,k
      real :: r,u(xp,tp),x,t,uy(10,10)      

      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,tp
         u( 1,j)=0.0
         u(21,j)=0.0
      enddo
      
      do j=1,tp-1     
         do i=2,xp-1     
            u(i,j+1)=r*(u(i+1,j)+u(i-1,j))+(1.0-2.0*r)*u(i,j)
         enddo
      enddo

      m=0
      l=0
      do j=1+k,tp,k
         l=0
         m=m+1
         do i=2,xp-10
            l=l+1
            uy(l,m)=u(i,j)
            t=(0.05/float(k))*float(j-1)
            x=0.05*float(i-1)
c        write(6,*) 'uy(',l,',',m,')', uy(l,m)
c            write(6,*) t,x,uy(l,m)
         enddo
c         write(6,*) ''
      enddo
                  
      return
      end      
