c  25. $B%]%"%=%sJ}Dx<0$rK~$?$9EE0L$N%0%i%U(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex25
c
      implicit none
      integer :: i,j
      real :: w,ep,pi,h,u(21,21)
c
c      write(6,*)''
c      write(6,*)''

      do i=1,21
         u(i, 1)=0.0
         u(i,21)=0.0
      enddo
      do j=1,21
         u( 1,j)=0.0
         u(21,j)=0.0
      enddo
c      do i=7,14
c         u(i, 1)= 1.0
c         u(i,21)=-1.0
c      enddo

      h =0.05
      pi=3.141592653589
      w =2.0-sqrt(2.0)*pi*sqrt((1.0/21.0)**2+(1.0/21.0)**2)
      ep=1.0e-3

      call SOR25(u,w,ep,h)

      stop
      end
c
c-----SOR$BK!$rMQ$$$?%]%"%=%sJ}Dx<0$r2r$/%5%V%k!<%A%s(B
      subroutine SOR25(u,w,ep,h)
      implicit none
      integer :: i,j,k
      real :: x,y,w,ep,Ee,h,f,u(21,21),up(21,21),Au,e(21,21)

      Ee=1.0
      do j=1,21
         do i=1,21
            up(i,j)=u(i,j)
         enddo
      enddo
   
      do while (Ee > ep)

       do j=2,20
          do i=2,20
             Au=(up(i-1,j)+u(i+1,j)+up(i,j-1)+u(i,j+1)-
     -           (h**2)*f(i,j))/4.0-u(i,j)
             up(i,j)=u(i,j)+w*Au
             e(i,j)=w*Au
          enddo
       enddo

       Ee=0.0 
       do j=2,20
          do i=2,20
             Ee=Ee+abs(e(i,j))
          enddo
       enddo

       do j=2,20
          do i=2,20
             u(i,j)=up(i,j)
          enddo
       enddo

      enddo

      do i=1,21
         x=0.05*float(i-1)
         do j=1,21
            y=0.05*float(j-1)
            write(10,*) y,x,u(i,j)
         enddo
         write(10,*)''
      enddo
      return
      end

c-----f(x,y)=sin(2pi*x)sin(pi*y)$B$N%U%!%s%/%7%g%sJ8(B
      function f(i,j)
      implicit none
      integer :: i,j
      real :: f,x,y,pi

      pi=3.141592653589
      x =0.05*float(i-1)
      y =0.05*float(j-1)
      f = sin(2.0*pi*x)*sin(pi*y)

      return
      end
