Hi John, If I understand correctly you are trying to visualize the square magnitude of some wavefunction of a 2D section of a 3D system. This problem has 2 parts: first the relevant wavefunctions components must be extracted from the wavefunction defined over the full (3D) system, and then these components must be plotted, e.g. by plotting a 2D colormap. The first part is relatively straighforward (and indeed you found a solution on the mailing list), however the second part has not been explicitly discussed on the mailing list and is indeed more challenging. First to tackle the mistakes that I saw in your script. In your system creation function you have the line return rsq < s^2 when what you probably want is return rsq < s**2 The caret (^) operator in Python means bitwise XOR. Normarily this would have raised an error, but as the type of 's' is 'int' in your example, and thus bitwise XOR is a well defined operation (though not the one you want). The second mistake is actually my fault, as the code on the mailing list that you linked to was actually written by me and contains a few typos. I realize now that we can simplify this code and make it more general: def in_sheet(site): return site.pos[1] == 0 # site in x-z plane, y position == 0 sheet = [idx for idx, site in enumerate(sysf.sites) if in_sheet(site)] # remember that 'wf(0)' returns a **sequence of wavefunctions** # so we first need to index which wavefunction we want (here the 1st one) wavefunction[0][sheet] Here we first construct a *list* of indices of all the sites in the x-z plane, and then index the wavefunction (a numpy array) with this list of indices. The second part, i.e. plotting the returned wavefunction, requires using some features of Kwant that, while documented, are rarely used directly by end-users. TL;DR the following code snippet will give you what you want: # remember that 'wf(0)' returns a **sequence of wavefunctions** # so we first need to index which wavefunction we want (here the 1st one) wavefunction_in_plane = wavefunction[0][sheet] density = np.abs(wavefunction_in_plane)**2 positions = [site.pos for site in sysf.sites if in_sheet(site)] # map positions to x-z plane positions_xz = np.array([(x, z) for x, y, z in positions]) data, min, max = kwant.plotter.mask_interpolate(positions_xy, density) image = pyplot.imshow(data, origin='lower') pyplot.colorbar(image) pyplot.show() We use the function 'mask_interpolate' from 'kwant.plotter' that will return data that can be directly shown as a colormap with 'pyplot.imshow'. We also need to map the site positions to the x-z plane, as 'imshow' will not know how to plot the output of 'mask_interpolate' if we provide positions in R^3 to this latter function. I have also attached a working example. Happy Kwanting, Joe