Problem linked here

I’m trying to plot the electric potential of a square area with finite charge density where the top and bottom edges have 0V and the side edges have 100 V. I’m asked to use the Gauss-Seidel method with overrelaxation (w in my code).

L=100
a=1
e0=8.85e-12
p=-1e6/e0
w=0.9
iterations=0
tolerance=1e-4
dvmax=100
vpot=np.zeros((L,L))
vpot[33:67,33:67]=p
vpot[:,L-1]=100
vpot[:,0]=100
while dvmax>tolerance:
      iterations+=1
      dvmax=0
      for i in range(1,L-1):
          for j in range(1,L-1):
              vnew=(vpot[i-1,j]+vpot[i+1,j]+vpot[i,j-1]+vpot[i,j+1]+(p*(a**2)/e0))/4
              dv=(1+w)*(vnew-vpot[i,j])
              if abs(dv)>dvmax:
                 dvmax=abs(dv)
              vpot[i,j]+=dv
print(f'Number of iterations: {iterations}')
plt.imshow(vpot)
plt.plasma()
plt.colorbar()
plt.show()
print(f'V(50,50) is {vpot[50,50]} V')

I’ve tried changing the sign of the charge density, the deltaVmax, and the spacing between points. I’ve also tried taking out the charge density but still including the poisson term in the code assuming that the code will calculate the charge density for me in the square area. Each variation runs for a very long time and never seems to end. I’m confused on what exactly is the problem.

I previously used essentially the same code where there’s no charge density, so p/rho is zero, and the plot came up quickly and without issue.