Program RELAX2
c
c  Sample program to solve Laplace equation using relaxation method
c
      Implicit NONE
      Integer XSIZE, YSIZE, XYSIZE
      Parameter (XSIZE=101, YSIZE=101, XYSIZE=XSIZE*YSIZE)
      Real MESH
      Parameter (MESH= 1.0/(XSIZE-1) )
      Real u(XSIZE,YSIZE) / XYSIZE*0.0 /
      Real EPS1,EPS2
      Parameter (EPS1=4., EPS2=2.0)
      Integer YBOUND
      Parameter (YBOUND=51) 

      Real corr
      Real cmax
      Integer LOOP
      Parameter (LOOP=30000)
      Integer n,i,j
c
c Initialization (ground potential for lower edgea, higher edge at 6 V)
c upper part EP2, lower part EP1
      Do i=1, XSIZE
         u(i,1)=0.0
         u(i,YSIZE)=6.0
         u(1,i)=0.0
         u(XSIZE,i)=0.0
      End do
c
c Relaxation
c
      Do n=1, LOOP
         cmax = 0.
         Do j=2, YSIZE-1
            Do i=2, XSIZE-1 ! Fortran array changes inner suffix first
c              corr=(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))/4 - u(i,j)
              if(j. eq. YBOUND) then
                corr=(EPS2*u(i,j+1)+EPS1*u(i,j-1))/(EPS1+EPS2)-u(i,j)
              Else
              corr=0.25*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-4.*u(i,j))
              Endif
              if(abs(corr) .gt. cmax) cmax = corr
              u(i,j) = u(i,j) + corr
c         floating condition for x edge
              if( i. eq. 2) u(1, j)=u(i,j)
              if( i. eq. XSIZE-1) u(XSIZE, j)=u(i,j)
            End do
         End do
c         if(mod(n,50).eq.0) Write(*,*) '# correction', cmax
      End do
c
c output phi array
c
      Write(*,*)'# READ MESH'
      Do i=1, XSIZE
          Do j=1,YSIZE
           Write(*,' (3f10.4)') real(i)*MESH, real(j)*MESH, u(i,j)
          End do
          write(*,*) ' '
      End do
      End