It's also straightforward to manipulate this algebra in python using the matrices I = np.array([[1,0],[0,1]]) #I=eye(2) X = np.array([[0,1],[1,0]]) Z = np.array([[1,0],[0,-1]]) Y = 1j*X.dot(Z) # Y = np.array([[0,-1j],[1j,0]]) H = (X+Z)/sqrt(2) x0 = np.array([1,0]) x1 = np.array([0,1]) Then use the kron ("kronecker direct product") function to define: def p(A,B,C,D,E): return np.kron(A,kron(B,kron(C,kron(D,E)))) and thereby operators on the 5 Qbit space M0 = p(I,Z,X,X,Z) M1 = p(Z,I,Z,X,X) M2 = p(X,Z,I,Z,X) M3 = p(X,X,Z,I,Z) I5 = p(I,I,I,I,I) H5 = p(H,H,H,H,H) X5 = p(X,X,X,X,X) and the encoded (0-bar and 1-bar) states b0 = (1/4.)*(I5+M0).dot(I5+M1).dot(I5+M2).dot(I5+M3).dot(p(x0,x0,x0,x0,x0)) b1 = (1/4.)*(I5+M0).dot(I5+M1).dot(I5+M2).dot(I5+M3).dot(p(x1,x1,x1,x1,x1)) Matrix elements are, e.g., b0.dot(H5).dot(b1) ================================================== Or in matlab, use the matrices I=[1 0; 0 1]; X=[0 1; 1 0]; Z=[1 0; 0 -1]; Y=1j*X*Z; H = [1 1; 1 -1] / sqrt(2) x0=[1 0]'; x1=[0 1]'; Then use the kron ("kronecker direct product") function to define: p = @(A,B,C,D,E) kron(A,kron(B,kron(C,kron(D,E)))); and thereby operators on the 5 Qbit space M0=p(I,Z,X,X,Z); M1=p(Z,I,Z,X,X); M2=p(X,Z,I,Z,X); M3=p(X,X,Z,I,Z); I5=p(I,I,I,I,I); H5=p(H,H,H,H,H); X5=p(X,X,X,X,X); and the encoded (0-bar and 1-bar) states b0=(1/4)*(I5+M0)*(I5+M1)*(I5+M2)*(I5+M3)*p(x0,x0,x0,x0,x0); b1=(1/4)*(I5+M0)*(I5+M1)*(I5+M2)*(I5+M3)*p(x1,x1,x1,x1,x1); Matrix elements are, e.g., b0'*H5*b1