Here’s what I eventually ended up with. Maybe you can see what you think. It has close to the same number of iterations by using the new suggested epsilon and putting fabs in the C code.
module rhofunc
implicit none
public
integer, parameter :: dp = kind(0.d0)
contains
pure real(dp) function rho(x, y)
real(dp), intent(in) :: x, y
if(x > 0.6_dp .and. x < 0.8_dp .and. y > 0.6_dp .and. y < 0.8) then
rho = 1.0_dp
else if(x > 0.2_dp .and. x < 0.4_dp .and. y > 0.2_dp .and. y < 0.4_dp) then
rho = -1.0_dp
else
rho = 0.0_dp
endif
end function
end module
program poisson
use rhofunc, only: rho
implicit none
integer, parameter :: dp = kind(0.d0), M = 100
integer :: i, j
real(dp), parameter :: epsilon0 = 3.85e-8_dp, target = 1e-6, a = 0.01
real(dp) :: delta, start, end, phiprime(0:M+1, 0:M+1), phi(0:M+1, 0:M+1), a2
integer :: iters
real(dp) :: rho_mat(M, M)
iters = 0
phi = 0.0_dp
delta = 1.0_dp
call cpu_time(start)
do i = 1, M
do j = 1, M
rho_mat(i, j) = rho(i*a, j*a)
enddo
enddo
a2 = a**2.0_dp
do while(.true.)
phiprime(1:M, 1:M) = (phi(2:M+1, 1:M) + phi(0:M-1, 1:M) + phi(1:M, 2:M+1) + phi(1:M, 0:M-1))/4.0_dp &
+ a2/4.0_dp/epsilon0*rho_mat
delta = maxval(abs(phiprime - phi))
iters = iters + 1
if(delta > target) then
phi = phiprime
else
exit
endif
enddo
call cpu_time(end)
print *, iters
print *, end - start
end program