!  A simple solution to the heat equation using arrays
!  and pointers

program heat2

real, dimension(10,10), target :: plate
real, dimension(8,8)           :: temp
real, pointer, dimension(:,:)  :: n, e, s, w, inside

real    :: diff
integer :: i,j, niter

! Set up initial conditions
plate = 0
plate(1:10,1) = 1.0  ! boundary values
plate(1,1:10) = (/ ( 0.1*j, j = 10, 1, -1 ) /)

!  Point to parts of the plate
inside => plate(2:9,2:9)
n => plate(1:8,2:9)
s => plate(3:10,2:9)
e => plate(2:9,1:8)
w => plate(2:9,3:10)

! Iterate
niter = 0
do
  temp = (n + e + s + w)/4.0
  diff = maxval(abs(temp-inside))
  niter = niter + 1
  inside = temp
  print *, niter, diff
  if (diff < 1.0e-4) then
    exit
  endif
end do

do i = 1,10
  print "(10f7.3)", plate(1:10,i)
enddo

end program heat2
