
      program main
      implicit none

c     Simple Jacobi relaxation for Laplace equation.
        
        integer n
        real eps
        parameter (n = 64, eps = 1.0e-6)

        real u (0 : n, 0 : n), uold (0 : n, 0 : n)
        real bdry
        integer i, j, niter
        logical convgd

c       Zero interior and initialise boundaries.

        do i = 0, n
          do j = 0, n
            if (i .eq. 0 .or. i .eq. n .or. j .eq. 0 .or. j .eq. n) then
	      u (i, j) = bdry (float (i) / n, float (j) / n)
            else
	      u (i, j) = 0.0
            end if
          end do
        end do

        convgd = .false.
        niter  = 0
        do while (.not. convgd)

c         Save old value.

          do i = 0, n
	    do j = 0, n
              uold (i, j) = u (i, j)
            end do
          end do

c         Update.

          do i = 1, n - 1
	    do j = 1, n - 1
              u (i, j) = (uold (i, j - 1) + uold (i, j + 1) +
     &                    uold (i - 1, j) + uold (i + 1, j)) / 4
            end do
          end do

c         Check convergence.

          convgd = .true.
          do i = 1, n - 1
	    do j = 1, n - 1
              convgd = convgd .and.
     &                 (abs (u (i, j) - uold (i, j)) .lt. eps)
            end do
          end do

          niter = niter + 1
        end do

        print *, niter
c        print *, u
      end

      real function bdry (x, y)
      implicit none
      real x, y

        bdry = x**2 - y**2
      end

