anchor.diff:  
   ! -- Boundaries and other initializations
   do i=1,n
  -  phi(i,m+1) = (real(i)-.5)*mesh_size
  -  x(i) = (real(i)-.5)*mesh_size
  -  a_n(i,m) = 2*a_n(i,m)
  -  a_s(i,1) = 2*a_s(i,1)
  +  phi(i,m+1) = (real(i)-.5_wp)*mesh_size
  +  x(i) = (real(i)-.5_wp)*mesh_size
  +  a_n(i,m) = 2.0_wp*a_n(i,m)
  +  a_s(i,1) = 2.0_wp*a_s(i,1)
   end do
   
   do j=1,m
  -  phi(n+1,j) = (real(j)-.5)*mesh_size
  -  y(j) = (real(j)-.5)*mesh_size
  -  a_e(n,j) = 2*a_e(n,j)
  -  a_w(1,j) = 2*a_w(1,j)
  +  phi(n+1,j) = (real(j)-.5_wp)*mesh_size
  +  y(j) = (real(j)-.5_wp)*mesh_size
  +  a_e(n,j) = 2.0_wp*a_e(n,j)
  +  a_w(1,j) = 2.0_wp*a_w(1,j)
   end do
   
  -phi(n+1,m+1) = 1.0
  -x(n+1) = 1.0
  -y(m+1) = 1.0
  +phi(n+1,m+1) = 1.0_wp
  +x(n+1) = 1.0_wp
  +y(m+1) = 1.0_wp
   
   write(axes_file,"(A,' ',A)")'x','y'
   do k=0,max(m+1,n+1)
  
positive.code: 
  ! ITERATE
  ! Goes through a maximum of 10000 iterations, but stops once
  ! RMS change is less than the specified value for three iterations.
  h=0
  write(convergence_file,"(3(A,' '),A)")'ITER','RMS','dRMS','H'
  do k=1,10000
    rms = 0.0
    do i=1,n
      do j=1,m
        delta = (relaxation/a_p(i,j))*(a_w(i,j)*phi(i-1,j) + a_e(i,j)*phi(i+1,j) + a_n(i,j)*phi(i,j+1) + a_s(i,j)*phi(i,j-1) &
          - a_p(i,j)*phi(i,j))
        phi(i,j) = phi(i,j) + delta
        e = abs(phi(i,j) - x(i)*y(j))
        rms = rms + e**2
      end do
    end do
    rms = sqrt(rms)
    if (abs(rms-rms_old) < error) then
      h = h + 1
    endif
    !write(convergence_file,"(I4,' ',2(F10.7,' '),I1)")k,rms,abs(rms-rms_old),h
    write(convergence_file,*)k,rms,abs(rms-rms_old),h
    if (h >= 3) then
      exit
    endif
    rms_old = rms
  end do
  
  ! WRITE SOLUTION
  
negative.code: call file_close(convergence_file)
  call file_close(phi_file)
  
  end program