Damn it, N is inverted and I noticed it now after posting. Sorry about that, here is correct one: ======================================== from numpy import arange, ones import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() ax = fig.add_subplot(111, projection='3d') o = ones(4) r = arange(4) # planes: for z in arange(3)+1: ax.bar(r, o*4, zs=z, zdir='x', alpha=.05, width=1) ax.bar(r, o*4, zs=z, zdir='y', alpha=.05, width=1) ax.bar(r, o*4, zs=z, zdir='z', alpha=.05, width=1) # N for i in [1, 2]: ax.bar3d([3-i], [0], [i], [.9], [.1], [.9], color='y', linewidth=.1) ax.bar3d(o+(i*(-1)**i), o-1, r, o-.1, o-.9, o-.1, color='y', linewidth=.1) # cage ax.bar3d([0], [0], [0], [4], [4], [4], alpha=.05, color='w', linewidth=0) plt.show() # plt.savefig('numpy.png') ========================================