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=2., EPS2=1.0)
      Integer YBOUND
      Parameter (YBOUND=81) 

      Real corr
      Real cmax
      Integer LOOP
      Parameter (LOOP=3000)
      Integer n,i,j
c
c Initialization (ground potential for all but 1 side)
c
      Do i=1, XSIZE
         u(i,1)=0.0
         u(i,YSIZE)=1.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
            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