Creating and updating 3 index matrix in Python

If I understand your question, you will need 2 changes:

1. First you need additional variables to check the positions that are fixed thru the iteration. I added sets with (j,k,l) tuples to do this. So you can follow my logic, I initially created 3 sets; 1 each for these indices: 1) fixed V1 (V1_set), 2) fixed V2 (V2_set) and 3) boundary (zero_set), then union all 3 sets into a single set (called check_set). You could start with a single set and update as you add. Side note: your code has V[1,:,:] = 0, but I think you really want V[0,:,:] = 0. Let me know if I interpreted that incorrectly.
2. Second, you need to loop on the axis length in each direction(attributes are V.shape[0], V.shape[1], V.shape[2]). Inside the loop I check each (i,j,k) against check_set, and only calculate anew V1[j, k, l] value if it is NOT in the set.

See code below:

V1 = 10
V2 = -5
Mx = 101
My = 151
Mz = 301

V1_set = { (46,k,l) for k in range(51,101,1) for l in range(101,201,1) }
V2_set = { (56,k,l) for k in range(51,101,1) for l in range(101,201,1) }

zero_set = set()
zero_set.update( { (0,k,l) for k in range(My) for l in range(Mz) } )
zero_set.update( { (100,k,l) for k in range(My) for l in range(Mz) } )
zero_set.update( { (j,0,l) for j in range(Mx) for l in range(Mz) } )
zero_set.update( { (j,150,l) for j in range(Mx) for l in range(Mz) } )
zero_set.update( { (j,k,0) for j in range(Mx) for k in range(My) } )
zero_set.update( { (j,k,300) for j in range(Mx) for k in range(My) } )

check_set = set().union(V1_set,V2_set,zero_set)

V = np.zeros((Mx, My, Mz)).astype(int)
V[46, 51:101, 101:201] = V1   #the values of these nodes should stay fixed throughout iteration
V[56, 51:101, 101:201] = V2   #the values of these nodes should stay fixed throughout iteration
V[1,:,:] =V[100,:,:] =V[:,1,:] =V[:,150,:] =V[:,:,1] =V[:,:,300] = 0     #the values of these nodes should stay fixed throughout iteration

for j in range(V.shape[0]):
    for k in range(V.shape[1]):
        for l in range(V.shape[2]):
            if (j,k,l) not in check_set:
                V[j, k, l] = (V[j+1, k, l] + V[j-1, k, l] + V[j, k+1, l] + V[j, k-1, l] + V[j, k, l+1] +V[j, k, l-1])/6

After posting the solution above, it occurred to me that the ranges used in zero_set are intended to avoid the first/last (array boundary) indices. If so, there is no need for zero_set. You can handle this by modifying the range arguments as shown below:

check_set = set().union(V1_set,V2_set)
for j in range(1,V.shape[0]-1):
    for k in range(1,V.shape[1]-1):
        for l in range(1,V.shape[2]-1):
            if (j,k,l) not in check_set:
                V[j, k, l] = (V[j+1, k, l] + V[j-1, k, l] + V[j, k+1, l] + V[j, k-1, l] + V[j, k, l+1] +V[j, k, l-1])/6

Additional observations to consider:

• I noticed you created array V with .astype(int). Are you sure that's what you want (and not floats)? In general, your calculations will not return integer values.
• The way your code is written, you are changing the values of V[j,k,l] as you go. So, you are using updated values of V[j,k,l] for j,k,l less than the current j,k,l, and previous V[j,k,l] values for j,k,l greater than the current j,k,l.
• Finally, I assume you are going to iterate thru this calculation until the change between 2 cycles is "acceptably small". If so, you need to have 2 copies of the array ("old" and "new") to take the difference. Take care to use .copy() when copying to create a new/different np.array object.

This is an updated answer based on new information and code added to initial post. You have at least 1 problem with your logic. The if (j,k,l) not in check_set: block skips over (j,k,l) values that you want to hold constant. As a result, you don't calculate Vnew at these points. That will cause problems calculation the change with each iteration (and will give the wrong result). Also, I think you need V = Vnew.copy(). Otherwise, V and Vnew reference the same object.

Here is my simple approach to iterate with a hardcoded error tolerance.

check_set = set().union(V1_set,V2_set)
Vi = V.copy()
Vn = np.zeros((Mx, My, Mz))
diff = max(abs(V1), abs(V2))
i = 1
print('Start Cycle#',i,'; diff =',diff)
while diff > 0.25:
    for j in range(1,V.shape[0]-1):
        for k in range(1,V.shape[1]-1):
            for l in range(1,V.shape[2]-1):
                if (j,k,l) in check_set:
                    Vn[j, k, l] = Vi[j, k, l]
                else:
                    Vn[j, k, l] = (Vi[j+1, k, l] + Vi[j-1, k, l] + Vi[j, k+1, l] + Vi[j, k-1, l] + Vi[j, k, l+1] +Vi[j, k, l-1])/6      
      
    diff = max(abs(np.amax(Vn-Vi)), abs(np.amin(Vn-Vi)))
    print('Cycle#',i,'completed; diff =',diff)
    i += 1
    Vi = Vn.copy()

This implementation will "converge" in 10 iterations. However, this only checks the error between two successive cycles is less than a hard coded tolerance (similar to the second part of the desired error check).
I did NOT implement the first error check: "use the maximum value of the local residual, defined as the (absolute value of the) difference between the potential value at the central node and the arithmetic average of the other values in the stencil." I am not 100 % sure of the intent. Is the stencil the 6 points around [j,k,l]? If so, I think you need a similar calculation AFTER you calculate the new Vn values, something like this:

error[j, k, l] = abs(Vn[j, k, l] - (Vn[j+1, k, l] + Vn[j-1, k, l] + Vn[j, k+1, l] + Vn[j, k-1, l] + Vn[j, k, l+1] +Vn[j, k, l-1])/6 )