What is the most efficient (quickest) way to multiply 20 identical 6x6 matrices (M)?

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = np.array([[w*p*q,w*q,0,0,0,0],
               [0,0,v,0,0,0],
               [0,0,0,nc,0,c],
               [0,0,0,0,v,0], 
               [w,w,v,nc,0,c],
               [0,0,0,n,0,1]])
for l in range(N-1):
    M = np.dot(M, Mi)
difZ = sy.diff(Z2,w)
expr = w*(np.divide(difZ,Z2))
Z_lamda = sy.lambdify([w,v,p,q], expr, "numpy")

For your special use case, I'd recommend using numpy.linalg.matrix_power (which wasn't mentioned in the linked question).

Timings

Here's the setup code I used:

import numpy as np
import sympy as sy
sy.init_printing(pretty_print=False)

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = M.copy()

and here's some timings comparing your original iterative dot approach to matrix_power:

%%timeit
M = Mi.copy()
for _ in range(N-1):
    M = np.dot(M, Mi)
# 527 ms ± 14.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%%timeit
np.linalg.matrix_power(Mi, N)
# 6.63 ms ± 96.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So matrix_power is about 80X faster.

Added bonus: matrix_power works better with arrays of Sympy expressions

For whatever reason, matrix_power seems to work better with Sympy than the iterative dot method. The expressions in the resulting array will be more simplified with fewer terms. Here's how you can count the terms in the result array:

import numpy as np
import sympy as sy

def countterms(arr):
    return np.sum([len(e.args) for e in arr.flat])

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = M.copy()

for _ in range(N-1):
    M = np.dot(M, Mi)

Mpow = np.linalg.matrix_power(Mi, N)

print("%d terms total in looped dot result\n" % countterms(M))
print("%d terms total in matrix_power result\n" % countterms(Mpow))

Output:

650 terms total in looped dot result

216 terms total in matrix_power result

In particular, print(Mpow) runs much, much faster than print(M).