how to fix python iteration for some 2D time dependent diffusion function using implicit schem

HI I'm new at this site and I'm looking for an answer for my python code

my code is simply solve groundwater head diffusion (using 2d implicit differentiation). in other words, I'm solving the function as time dependent
with Laplace Equation : d^2h/dx^2+d^2h/dy^2 = ds/dt where dh is the change of head over flow path (dx = dy) in space ds is some constant value (storage coefficient) and dt is the time increment

so I do solved the example under the condition of to boundary conditions upstream 5 m and down stream of 3 m and my grid is 7x7 cells

I solved the problem using C# and matlab and both answers were close after conversions ( which reach the equilibrium for some exact analytical solution) but when I used Python my solution didn't reach the equilibrium

here is my code with python

 import numpy as np

h = np.zeros((7,7,10))
h0 = 5.
K = 2.
DT = 1.
DX = 100.
S = 0.0005
B = 8.
D = (K*B)/S
r = (D*DT)/(DX*DX)

dummy = h.shape
nrow = dummy[0]
ncol = dummy[1]
ntm = dummy[2]
print 'Groundwater Head matrix is a ', nrow, ' by ', ncol, ' matrix.','with',ntm,'increments'
for i in range(nrow):
    for j in  range(ncol):
        for k in range(1):
         h[i,j,k] = h0
         continue

        h[0:,0,1:] = 5.
        h[0:,6,1:] = 3.

        ni = 1
conv_crit = 1e-9
converged = False
while (not converged):
  max_err = 0
  for i in range(1, nrow - 1):
      for j in range(1,ncol - 1):
          for k in range(1,ntm-1):

            h_old = h[i,j,k+1]
            h[i,j,k+1] = (r*h[i-1,j,k+1]+r*h[i+1,j,K+1]+r*h[i,j-1,K+1]+
                          r*h[i,j+1,k+1]+h[i,j,k])/(1+4*r)

            h[0,0:6,k+1] = h[2,0:6,k+1]
            h[6,0:6,k+1] = h[4,0:6,k+1]

            diff = h[i,j,k+1] - h_old
            if (diff > max_err):
             max_err = diff
  if (max_err <= conv_crit):
      converged = True
  ni = ni + 1
  print 'Number of iterations = ', ni - 1
  print h[0:,0:,9]
  print ' Error =', max_err

this is my result with python

Number of iterations =  132
[[ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.     ]
 [ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.        ]
 [ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.        ]
 [ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.        ]
 [ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.        ]
 [ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.        ]
 [ 5.          4.40796148  3.88890558  3.5072111   3.24606331  3.0865587  3.        ]]
 Error = 9.17176556925e-10

and here with matlab

5.0000    4.6599    4.3229    3.9892    3.6583    3.3291    3.0000
    5.0000    4.6592    4.3217    3.9879    3.6573    3.3285    3.0000
    5.0000    4.6599    4.3229    3.9892    3.6583    3.3291    3.0000
    5.0000    4.6606    4.3240    3.9904    3.6593    3.3296    3.0000
    5.0000    4.6613    4.3250    3.9915    3.6602    3.3301    3.0000
    5.0000    4.6619    4.3259    3.9925    3.6610    3.3305    3.0000
    5.0000    4.6613    4.3250    3.9915    3.6602    3.3301    3.0000

matlab use only 41 iteration

Notice that I used the Gauss Seidel Iteration