I am trying to plot the wavefunction for a 2D NISIN Hamiltonian however I am unable to get a result, this is my code: import kwant import tinyarray as ta import numpy as np import matplotlib.pyplot as plt from tqdm import tqdm import winsound import os from datetime import datetime import scipy.sparse.linalg as sla sig_0 = np.identity(4) s_0 = np.identity(2) s_z = np.array([[1., 0.], [0., -1.]]) s_x = np.array([[0., 1.], [1., 0.]]) s_y = np.array([[0., -1j], [1j, 0.]]) tau_z = ta.array(np.kron(s_z, s_0)) tau_x = ta.array(np.kron(s_x, s_0)) tau_y = ta.array(np.kron(s_y, s_0)) sig_z = ta.array(np.kron(s_0, s_z)) tau_zsig_x = ta.array(np.kron(s_z, s_x)) tau_zsig_y = ta.array(np.kron(s_z, s_y)) # Lorentzian definition def lorentzianxy(x, y, ind,y0): gam = 3 L = params["L"] #if 0L/2: return (4 * t - Smu) * tau_z + B * sig_z + Delta * tau_x - potential(site, pot,ind,y0)*tau_z # else: #return (4 * t ) * tau_z + B * sig_z + Delta * tau_x - potential(site, pot,ind,y0)*tau_z def onsiteNormal(site, pot,ind,y0): return (4 * t - mu) * tau_z #- potential(site, pot,ind,y0)*tau_z def onsiteBarrier(site, pot,ind,y0): return (4 * t - mu + barrier) * tau_z #- potential(site, pot,ind,y0)*tau_z def hop_x(site, pot, ind,y0): return -t * tau_z + 0.5j * alpha * tau_zsig_y def hop_y(site, pot, ind,y0): return -t * tau_z - .5j * alpha * tau_zsig_x # On each site, electron and hole orbitals. lat = kwant.lattice.square(norbs=4) syst = kwant.Builder() # S syst[(lat(i, j) for i in range(1,L-1) for j in range(W))] = onsiteSc barrierLen = 1; # I's syst[(lat(i, j) for i in range(barrierLen) for j in range(W))] = onsiteBarrier syst[(lat(i, j) for i in range(L-barrierLen, L)for j in range(W))] = onsiteBarrier syst[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x syst[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y # N's lead = kwant.Builder(kwant.TranslationalSymmetry((-1,0)), conservation_law=-tau_z#, particle_hole=tau_y ) lead[(lat(0, j) for j in range(W))] = onsiteNormal lead[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x lead[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y syst.attach_lead(lead) syst.attach_lead(lead.reversed()) return syst def sorted_eigs(ev): evals, evecs = ev evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose())))) return evals, evecs.transpose() params = dict(mu=0.3, Delta=.1, alpha=.8, t=1.0, barrier = 2.0, pot = 0.0,W = 5, L = 30, ind = 0, B = 0.3, Smu=0.00,y0=0) sys = makeNISIN2D(params) sys = sys.finalized() # Calculate the wave functions in the system. ham_mat = sys.hamiltonian_submatrix(sparse=True, params = params) evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0)) # Plot the probability density of the 10th eigenmode. kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1) And the error I get is this: "ValueError: The number of sites doesn't match the number of provided values." When I print the values they indeed are more than what they should be. I've tried a couple of things but I'm quite at a loss of what to do. I want to plot the wavefunction so I can see the Majoranas at the edges of the superconductor so that I can see how they respond to a local perturbation in the potential. Any tips would be very helpful. Kind regards, Hanifa