Program RELAX1
c
c  Sample program to solve Laplace equation using relaxation method
c
      Implicit NONE
      Integer XSIZE, YSIZE, XYSIZE
      Parameter (XSIZE=51, YSIZE=51, XYSIZE=XSIZE*YSIZE)
      Real MESH
      Parameter (MESH= 1.0/(XSIZE-1) )
      Real u(XSIZE,YSIZE) / XYSIZE*0.0 / 

      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
              corr=(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))/4 - u(i,j)
              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