!  A simple solution to the heat equation using arrays
!  and section subscripts

program heat1

real, dimension(10,10) :: plate, temp

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 ) /)

! Iterate
niter = 0
do
  temp(2:9,2:9) = (plate(1:8,2:9) + plate(3:10,2:9)        &
                  +plate(2:9,1:8) + plate(2:9,3:10))/4.0

  diff = maxval(abs(temp(2:9,2:9) - plate(2:9,2:9)))
  niter = niter + 1
  plate(2:9,2:9) = temp(2:9,2:9)
  ! To show how the convergence is progressing
  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 heat1
